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ABSTRACT 



Integral equations (IE's) are widely utilized to calculate induced currents on anten- 
nas and scatterers, but they are seriously restricted in their ability to handle inhomoge- 
neous penetrable structures having multiwavelength dimensions. The utilization of finite 
element (FE) techniques has not been as pervasive as the use of IE's. The IE represen- 
tation matrix is "full", containing few, if any, zero valued elements. The techniques for 
operating on these large-sized full matrices require undesirable amounts of processor 
time. FE techniques produce sparse matrices due to the strictly local interactions be- 
tween discrete unknowns. The application of FE's to unbounded problems, however, 
requires supplementarv’ enforcement of the far-field radiation conditions. The Field 
Feedback Formulation (F^) circumvents the full-matrix computational "bottleneck" by 
allowing FE based numerical methods to be employed. Even though the resultant sparse 
matrices may be larger than the "full" matrices discussed earlier, most elements have a 
value of zero. Numerical procedures exist to optimize operations with these sparse ma- 
trices. Calculational speeds can be orders of magnitude faster. Computer techniques to 
implement and validate this new technique are the basis for this thesis. Excellent 
agreement with classical results are demonstrated. 
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I. INTRODUCTION 



A. HISTORY 

The application of finite element techniques to evaluate the solution of differential 
equations is well documented. The utilization of these techniques by the 
electromagnetics community has not been as pervasive as the use of integral equations 
(IE's). IE's are widely utilized to calculate induced currents on antennas and scatterers. 
but they are seriously restricted in their ability to handle inhomogeneous penetrable 
structures having multiwavelength dimensions. As the complexity and number of nodal 
degrees of freedom grow, the size and dimension of the representative matrix must also 
grow. This matrix is "full", containing few, if any, zero-valued elements. The available 
numerical techniques for operating on these large-sized full matrices require undesirable 
amounts of processor time. 

Differential equation (DE) based techniques, such as the finite element method, 
produce sparse matrices due to the strictly local interactions between discrete unknowns 
which result. The application of DE's to unbounded problems, such as those of scat- 
tering and radiation, require some form of supplementary enforcement of the proper 
far-field conditions. These radiation boundary conditions are innately incorporated into 
integral equations. 

B. FIELD FEEDBACK FORMULATION 

The Field Feedback Formulation (F^) circumvents the full-matrix "bottleneck" in the 
computational process by allowing DE based numerical methods to be employed. Even 
though the resultant sparse matrices may be larger than the "full" matrices discussed 
earlier, most elements have a value of zero. Numerical procedures exist to optimize 
operations with these sparse matrices. Calculational speeds can be orders of magnitude 
faster than with full matrices (Ref 1]. Although the employs sparse matrices to rep- 
resent the fields in the materials being considered, it does require augmentation to en- 
force the radiation condition at infinity on the scattered fields. This comes in the form 
of a feedback matrix composed of surface integration generated elements. A concept 
evaluation, for a special axisymmetric case, was already accomplished, as detailed in 
[Ref 2] and (Ref 3]. Computer techniques to implement and validate this new technique 
are the basis for this thesis. 
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C. POTENTIAL BENEFITS 

This thesis will lead to an increased understanding of the advantages and disadvan- 
tages of this novel computational procedure for handling geometrically complex material 
scatterers. This method may ultimately allow computer-aided design of important 
electromagnetic structures such as low-observable aircraft, high efficiency dielectric lens 
antennas, and other electromagnetic scattering occurrences due to atmospheric anoma- 
lies. Structural details and material inhomogeneities, as well as physical dimensions (in 
multiple wavelengths), can be accommodated using the Field Feedback Formulation. 
These capabilities far surpass those which are possible with contemporary integral 
equation techniques for the case of inhomogeneous penetrable scatterers and antennas. 






II. FORMULATION 



A. INITIAL NOMENCLATURE 

Assume there is a three dimensional object that is infinite in one direction. Such an 
object, in cross section could look like Figure 1. This object only varies in two dimen- 
sions, and therefore, is actually a two dimensional (2-D) object. 

y 



Figure 1. 




The wavenumber is defined as 



, _ 2n _ 

*»- - c • 

where is the free space wavelength associated with an electromagnetic wave of fre- 
quency fg and c is the speed of light. The X and Y Cartesian coordinates are 



3 



wavenumber normalized, such that X = and k,^’. This coordinate normalization 
will be used throughout this development. A similar technique could also be developed 
in the polar coordinate sys tem . The magnetic field will also be normalized such that 
II = where i/o = ■sj'^ = I20n « 311LX is the impedance of free space and 

is the usual magnetic field in units of A/m. Thus the normalized 7/ has the same V/m 
units as E. Potentials may be defined as, 

E,{x^) = il,,{X,Y) ( 2 . 1 ) 

for the transverse magnetic (TM) case, with Ey, and 7/^ = 0, and 

= ( 2 . 2 ) 

for the transverse electric (TE) case, with 77„ Hy and E^ = 0 . 

Our objective is to calculate the scattered fields for an arbitrary (2-D) penetrable 
object using the Field Feedback Formulation (F^). As shown in Figure 2, a familiar 
closed loop system illustrates the relationship between the incident, scattered, total and 
far fields. 




Figure 2. Field Feedback Foriiiulation 

At point 1, the incident field drives the total system. I he incident field at 1, combined 
with the scattered field at 4, forms the total field at 2. This total field forms the 
boundary conditions that drive the finite element program at 2. At point 2, initially, the 

. 3 

incident field drives the U operator. U represents the feed forward operator in the F . 
This operator uses a finite element technique to solve the boundary value problem. At 
point 3, the boundary conditions are solved for the object perimeter potentials and the 
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normal derivative of these potentials. T represents the field feedback operator that takes 
the perimeter potentials and associated derivatives and provides the scattered fields at 
point 4, the offset boundary. The fields at point 1 and 4 are added (unlike the familiar 
feedback or control system where negative feedback is employed). At point 2, a com- 
bined incident and scattered field exists. These fields are the combined boundary con- 
ditions for the finite element boundary' value program. These combined fields (total 
fields) are then applied to the U operator to calculate the perimeter fields and the asso- 
ciated derivatives on the objects perimeter. This looping may be repeated until a steady 
state condition, at point 3, is reached. The existence of a steady state condition assumes 
stability. Stability for physical systems should not be a problem. However, when 
mathematically modeled, instabilities may result. The error magnification or condition 
number of the system must also be seriously considered if this iterative looping process 
is to be used. The alternative approach is to form an equivalent system, where: 

equivalent operator = [/ — J. 



with 



I — Identity Matrix 



and 

T • U = Combined Effects of the T and U operators . 

Either approach is viable since. 

^ total ^ incident 7"» L • ^ incident (E* L") • incident ... — [.! 7 • ^incident 

This becomes more obvious when is factored from the equality leaving, 

\ + T*U+iT*lf + . . . 

Placing this in closed form, 

OO 

^(r. IT'- 

n=0 
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Finding the equivalent operator will require a matrix inversion and for ver\' large prob- 
lems this may lead to excessive computation times. The matrix inversion technique will 
be investigated. The fields may then be extended to any point in space using a far field 
Green's function surface contour integral. This far field pattern is available at point 5. 

The incident field is usually produced by a plane wave generator. This provides the 
boundaiy conditions on the offset contour. This contour is called "offset" since it is 
approximately the same "shape" as the objects perimeter but is slightly larger. The dis- 
tance between the perimeter and this contour is called the offset distance and will be 
discussed in Chapter III. The boundaiy conditions may be any desired field or wave that 
satisfies Maxwell's equations. These waves may arrive from any direction and be of any 
magnitude. The boundary conditions may also be a composite of any number of waves 
since superposition does apply to these systems. A user provided subroutine is necessaiy 
if conditions other than a single plane wave, cylindrical mode {of arbitraiy mode num- 
ber) or individual input boundaiy condition is desired. 

B. MAXWELL'S EQUATIONS 

Maxwell's equations can be written using our previous normalizations as, 

VxE = HrH (2.3) 

and 

V X 77 = c,E. (2.4) 

[Ref 4 ] Let . with similar definitions for Dy and Dy. Equations 2.3 and 2.4 

can be further expanded into the Dy, D^and Dy components such that. 



IJ-yHy. = DyE. 


(2.5) 


HfHy = —DxEy 


(2.6) 


IJ-yHy = D^Ey DyEj. 


(2.7) 


= Dylly 


(2.8) 


tfEy = —DxIIy 


(2.9) 


CfE, = DxHy —DyH^. 


(2.10) 



6 



For the TM case, with propagation in the z direction, 






Mr 



and 



//,= 



-DxE, 

Mr 



Note that the H, field = 0. These two equations can be combined to form. 




( 2 . 11 ) 



Similarly for the TE case, with propagation in the z direction. 






Substituting equations 2.5 and 2.6 into equation 2.10 yields, 



. -E>xE, 



Substituting equation 2.1 into equation 2.13 yields. 



«r'Ai + + ^>( '^1 



= 0 . 



( 2 . 12 ) 



(2.13) 



(2.14) 



Equation 2.14 can be further simplified to, 

V* + =0. (2.15) 

Similarly, equations 2.2, 2.8, 2.9 may be substituted into equation 2.7. This yields, 

^•[-t^'A2]+Mr'A2 = 0. (2.16) 

Equations 2.15 and 2.16 are TM and TE duals. These two differential equations describe 
the potentials inside the object of interest. Defining = a and £, = /? for the T.M case 
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and -^ = a and pl, = fi for the TE case and substituting these new definitions into 
equations 2.15 and 2.16 yields one differential equation, 

V • = 0. (2.17) 

C. VARIATIONAL EQUIVALENCE TO THE DIFFERENTIAL EQUATION 

The Euler-Lagrange variational formulation is based on the stationarity of a func- 
tional. of the function and its first derivatives. [Ref 4 ] 



-JJ 

inside S 



F{X, Y, <A, Vi//) dXdY 



(2.18) 



It can be shown that the first variation of the functional is zero, <57 = 0 , if the 
Lagrangian. F, satisfies the Euler - Lagrange equation; 



c 

ex 



dF 



a(Av'A) 



BY \ d{DyiP) 




(2.19) 



The problem thus becomes to find the F, which when substituted into equation 2.19, 
yields the original differential equation. The Lagrangian, 

F= a{(7)^-.A)' + (7)>..A)'} - /?iA' 



can be simplified to. 



F= aVi// • Vi// — /?i//^. 



( 2 . 20 ) 



When equation 2.20 is substituted into equation 2.19, 



BF 



= 2a.D^->p 



cF 

BiDyiP) 



2aDy\p 



8F_ 

Cl// 



= 2IBrp. 
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Therefore, 



-I- {2uD,-^P) + 4^ {2y.Dy^) + 2/?<A = 0 

C*^ 1 O I 



or, 



Z);^-[aZ)_Y<^] + Dy[aDylp] + pip =0 



which simplifies lo, 



V . [aViA] + = 0. (2.21) 

Therefore, the general functional has been found since equation 2.21 and equation 2.17 
are identical. Substituting the a and /? definitions into equation 2.20 yields, 

( 2 - 22 ) 

and 

= (2.23) 

Equations 2.22 and 2.23 are integrated over the interior to S with known boundary 
conditions for either \j/^ or \p 2 on S. To physically interpret the variational formulation, 
it is noted that, 






and 

V,A2 = £,{£;,}"- £yAj- 

Therefore, 



and 



A 



H,H • H -c,E • E dXdY 
s 
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€,E • E - • H ciX dY. 



Substituting Maxwell's equations into /, gives, 






(V X £).//-(V X H)*EdXdY 



%) *) 

s 






V.(£ X H)dXdY 
s 



h = 



E X H^ ndl. 



as 



Note that £ x // is the complex oscillatory Poynting vector. It is dilTerent from the 
usual Poynting vector of £ x H'. Thus, both functionals, /, and /j are proportional to 
the complex phasors for oscillatory power. The oscillatory power is the phasor repres- 
enting the excess of instantaneous radiated power minus the average radiated power. 

D. FINITE ELEMENT BOUNDARY VALUE SOLUTION 

With the discussion of the variational mathematics completed, a simple rectangular 
boundary value problem will be discussed in detail. The rectangular geometry allows for 
an easier formulation but in no way limits the solution from being extended to more 
complicated object geometries. Figure 3 on page 11 shows the region of concern. 
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b 



>•)! 

CaU, y)) 



in S 



Boundary 

dO. 



□ ' ' -»• 

0 q X 

Figure 3. Rectangular Object 

Consider spanning the rectangular region by a triangular mesh, as shown in 
Figure 4 on page 12. Both a and /? may be functions of position but may not vary 
within an individual triangular element. Given a fine enough mesh structure, a smooth 
transition in material properties may be approximated. For notational simplicity this 
positional dcpendance will not be carried for\^ard. 1 he variational approach will yield 
the solution to the boundary value problem by finding the rp{x,y) which gives the sta- 
tionary value of. 



/ = 



ra 



(aViA • ViA - dx dy 



»'0 »'0 



which is constrained on the boundary by the previously specified boundary conditions. 
[Ref 5 ] 



N+1 




1 he values of ij/ at the interior nodes become the discretized unknowns: 

\l/jj = ^{Xf, Yj) /= 1. . . /)/ and j = \ ... N 



where 



Xi = i»AX 



and 



Yj=J.AY 



for the uniform mesh structure with AX = 



A/+1 A'+1 



and Ar = 



(/V+1) 



. Approximating, 






1=0 7=0 



which includes the known boundary nodal values \]/^j and for j = 0 and N+1, 
and i/^, 0 and for i = 0 and M + 1, linear pyramidal basis functions, u,j{x,y), which 
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have unit value at the (i,j) node and zero value at all surrounding nodes. Figure 5(a) is 
a top view of the pyramidal basis function, while Figure 5(b) is a perspective view. 





Figure 5. Pyramidal Basis Function, (a) top view, (b) perspective view 



1 he functional I will thus be a discrete function of each of the nodal values of 

^ = ^('Ao.O' •Ao, 1» •••' 'A/J* •••' 'Aa/+ 1, Af+l)- 
The approximate discrete solution will be found by the system, 

■ — = 0, for m = 1 . . . M and m = 1 ... A'. 



Now, 



ei 




8V<j/ 






dij) 

^^m.n 




dx dy = 0 
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where. 



A/+1 A'+l 



V 'A (-V, >■) = ^ ^ 'A,;V Uiji-X, >•) . 
i=0 y=0 



The gradient of the basis function is, 






dVrJ/ 

n 



and the basis function is, 






34 / 






Therefore, the system of equations to solve becomes. 



*^0 



V 2 /„, „ • V 1/2 — /? 2 /„, „\j/) J.x dy = 0. for m= \ . . . M and « = 1 . . . N. 



Substituting for Vi/' and 1/2 in terms of the nodal values of 1 // for 
2/1=1... M and « = 1 . . . A’ sives. 



A/ 4-1 V4-1 



jzi Jo Jo 



where. 




Vuij - „Ui j) dx dy = F{{tn, n),{ij)]. 



Regrouping this to put the known nodal values on the right hand side gives for 
/7z = 1 . . . .1/ and n = 1 . . . A' (interior nodes). 



M S 

/=! y=i Boundary 

Sodes Only 
for (if) 
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By renumbering the nodes using a single index for the unknown nodal values, 
k = S{i — 1) +7 and / = S{m — 1) + az and, 

M-S 

k=\ k' 

The functional F{1, k) = 0 if nodes / and k are not both associated with at least one 
common triangular element. Therefore, Vn, • Vaa* and u,u^ will be zero except in triangles 
where / and k both appear as nodes. For the example mesh structure of Figure 3 on 
page 11, this produces a banded matrix, F. This sparse matrix can be easily displayed 
by first defining column vectors of the unknown nodal values in the mesh 

The F - matrix elements F [(/ai, A7).(/,y)3ire zero unless the (m.n) node shares at least one 
element with the (i.j) node. Thus, the node values in \p, will be coupled only to 
<A,_i and any associated boundary nodes. This is written as. 



where A^, 5, and C, are .V x A' complex arrays and P/ is a A' x A'^^ array where A'js_ is the 
number of boundar\' nodes associated with the i-th column. This appears as. 



- 








- 


- 




- 


~ ~ 




c, 


















B2 








^2 




Pi 






^3 


^3 C3 


0 




• 




Pi 

• 


• 




0 








• 




• 

P\l-\ 


• 










Bm 


^A/ 




Pxf 




- 








- 


_ 




- 


_ _ 



If we denote i/aq as the initial boundar}’ values and as the final boundary^ values 
then \{/^^ becomes only the boundary conditions on the top and bottom at J ^0 and 



15 



iV+ 1. Note that the above system is tri-block in nature and has a large number of zero 
valued elements. The zero valued elements were omitted for clarity. Bach element is 
actually a matrix and therefore it is evident how sparse this system is. Each of the 
/1„ Z?,. Cl and P, matrices are equally sparse, however, a global synunetry does not ap- 
pear. 



E. EVALUATION OF THE F - MATRIX CONTRIBUTIONS 

Given an arbitrary element as shown in Figure 6, the potential, ip, can be linearly 
approximated by (Ref 5 ], 



<A(-v,>-) = (.t,r,i).[r] 




<A3 



1 




Figure 6. Mesh Element 



where i//* = is the nodal value at the k-th node, [T] = 3x3 Transform Array 

and 
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Ui.{x,y) = Linear basis functions for the k-th node 



= (^. 1) 



Tuk 

Ti,k 



= Ti,k- 



Note that, 



^k^Xm^ym) 



1 , for k = m 
0 , for k i=- in 



It can be shown that. 
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where \A\= triangle area, and 



2/1 = det 



X2 >2 



•^3 >‘3 



1 

1 

1 



2A = (x;^V3 + X3^V, + XiV-2) - (X3^V’2 + ^i>h + .v^yi) 



Furthermore, the linear basis functions, u*(x,j')> can be interpreted as relative areas of 
the triangle shown in Figure 7 on page 18. 
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Figure 7. Basis Functions Interpretation 

/4 = /I] + A 2 + 

is constant as varies and, 



Within a given triangular element, the evaluation for k = 1, 2, 3 and / = 1, 2, 3 in the 
element is of interest and, 



The assumption is made that a and /? will be approximated as constants within the tri- 
angle. These material constants can, however, vary from element to element. Taking 
the gradient terms first. 




inside 

triangle 

q 



^^k=T\.kX + T2^ky 




18 



Therefore. 



J J a^Vui^Vnndxdy 

triangle 



where is the area of the q-'‘ element. Next it can be derived that, 



r 



UiUi^dxdy = 



triangle 



f^MI 

H'-'i 



for k ¥= t 
fork = I, 



More generally, with each u, raised to an integer power. n„ 



c c 

ul'iifu^^dxdy = 2\A\ 

%} 

triangle 



>hW>h'- 

( a 7 j + th + + 2 )! 



The final result is, 



Fqih = J 1 - !^uiuf}dxdy 

triangle 

q 

= k = L 



F. GREEN'S FUNCTION CONTOUR INTEGRAL 

The scattered fields, \J/, from an arbitrary object in a vacuum satisfying Helmholtz's 
equation (see Figure 8 on page 20 ), 

vV + ^V = o. 



are [Ref 6 ], 
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Figure 8. Green's Function Integration 



where the Green's function is, 



and 



dn 



= n.Vil/ 



on the contour 



and 



dG 



on 



A A J^O 

= n • r — r~ 



H?{ko\r-F\). 
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The Hankel functions of the second kind of order zero and one, and , will 

present a problem for the numeric integration discussed in Chapter V. The imaginary 

portion of these functions rapidly approaches negative infinity as the argument ap- 

d\!/ 

proaches zero. The — is obtained bv a finite difference method usins the field 

cn ' 

boundary conditions on surface B (boundary conditions) and the calculated field condi- 
tions on surface P (object perimeter). This results in, 



Olj/ V boundary Y perimeter 

cn offset distance 



It will be shown that to maximize the accuracy of the Green's function the offset 

ciA 

distance should be made as large as possible. This, however, causes the — to be in- 



cn 



accurate. Thus, an optimal condition must be found that maxintizes the accuracy of the 
entire numeric integration. Such a condition does not maximize the accuracy of any one 
of the contributing parts to the Green's function integrand. 



G. FAR-FIELD EVALUATION 

When the Green's function integral discussed above is used for far field calculations 
several simplifying relationships develop. These simplifications require a less demanding 
numerical integration. To be in the far field region three conditions must exist, 

I F 1 > Z). I F I > Zq and I F | > , 

''■0 



where Zp is the free space wavelength and D is the maximum dimension of the object. 
[Ref 7 ] As X approaches infinity. 
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-Jx 



and 







This requires, 



C(f|p)=| 







and 
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In the far field 




cn A 



Tzk^R 



\! -y and rt • y? -► rt • f. Thus, 

g-J>v . 



Therefore, 



G(F|?')=^ 



2; 






TrAr/ 



and 



"'''(Flr) = -4 



C// 



0 * 

COS ^ 

-T— ^ (// • r)e ° 



nk^r 



With these new definitions substituted into the original Green's function integral 
equation. 



•A, 



scaaerecv 



lr) = 



J 



^nky 



-jk^ 



J -J- + Aq/2 •r\]/{r) 

C /7 



^'Aor' cos ^ 






Note that this equation is partitioned into a distance dependent term and a theta depend 
term. The theta dependent term may be defined as 






/ = 



Jr 



dip 



on 



+ k^n • rp{r 






COS 0 



dc’. 



The two dimensional bistatic radar cross section (RCS) per unit length of the cylindrical 
structure may now be defined as. 



RCS = a((p\ 4>‘] = lim 

^ / r-oo p‘>ic 



= lim 

r—^oo 



2nr\f\^ 

ui' 



22 



where the wavenumber, 
magnitude, \ij/‘\ = 1.0. 



a = lim 2nr 






I/I' 



Snhr 4A'n 



An = ■ 



and the incident field is assumed to be of unit 
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III. MESH GENERATION 



A. INTRODUCTION 

A display or plot of the computer generated mesh structure is not required for the 
problem solution, however, it provides an immediate visual confirmation that the in- 
tended problem geometry’ has been entered correctly into the computer. A large amount 
of initialization data is required for even the most rudimentary problem. For this reason, 
the input data is provided to the mesh generation program via a data file. Modifications 
are possible at a later time with minimum ellbrt. 

B. INPUT DATA 

The input data file is called INPUT.DAT and contains 25 input fields. Of these 
fields, 14 relate directly to the generation or display of the mesh structure. Only the 
object surface coordinates require adherence to a specific format. All other data need 
only be of the correct type (i.e., character, integer or real). A brief description of each 
field is provided below. 

• Field I is a label of no more than 12 characters. This label is for the plot and input 
data file. 

• Field 2 is a character flag that specifies the input coordinate system. If set to "R" 
or "r", the rectangular system is used. If set to "P" or "p", the polar system is used. 
If a "P"."p" ."R" or "r" is not detected, an error is returned to the display. 

• Field 3 is a character flag that if set to "I" or "i" will cause several intermediate 
values to be stored to disk during the Finite Element Boundary Value (FEBV) 
program execution. This option was used during the debug process. 

• Field 4 is a character flag that if set to "D" or "d" will create a DISSPLA 
FORTIUAN program capable of replicating the input object, in wavenumber nor- 
malized coordinates, on mainframe computers having a DISSPLA graphics pack- 
age. DISSPLA is a subroutine-based language. The generated program, called 
DISSPLA. FOR, is a compilation of four subroutines calls per element. This file 
can get very large for dense mesh structures. 

• Field 5 is a character flag that if set to "U" or "u" will cause a uniform material to 
be assumed. In the uniform case, no material interface exists. This option was 
only used to verify the Finite Element solution accuracy. 

• Field 6 is a character flag that if set to "M" or "m" will cause only the mesh to be 
generated. This is very useful when first starting a problem and the optimal mesh 
structure has not been determined. 

• Field 7 is a real number specifying the desired mesh resolution in wavelengths. The 
mesh resolution determines the dimension of the mesh elements. 
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• Field 8 is a real number specifying the distance, in wavelengths, between the object 
perimeter and the offset boundary contour. 

• Field 9 is a real number that specifies a multiplicative scaling factor for the nu- 
merical integration stepping function discussed in Chapter V. 

• Field 10 is a bias term that can be used to shift the numerical integration stepping 
function. This term is used for distances less than 1.0. 

• Field 1 1 is a bias term that can be used to shift the numerical integration stepping 
function. This term is used for distances greater than 1.0. 

• Field 12 is a real number specifying the maximum distance beyond which no further 
contribution to the Green's Function Integral is made. If this feature is not de- 
sired. this term should be made larger than the objects maximum dimension plus 
twice the offset distance. 

• Field 13 is an integer specifying the number of input data points. 

• Field 14 is an integer specifying the angular resolution, in degrees, desired for the 
final radar cross section calculation. 

• Field 15 is an integer specifying the mesh generation technique. 

• Field 16 is an integer specifying the perimeter node from which the bisection seg- 
ment originates. This node is called the "start node". 

• Field 17 is an integer specifying the perimeter node on which the bisection segment 
terminates. This node is called the "stop node". 

• Field IS is a pair of real numbers (on two lines) specifying the x and y coordinates 
by which the object will be displaced. 

• Field 19 is a real number specifying, in wavelengths, the desired distance between 
the origin and the first input data point. Fields IS and 19 when used together, al- 
low an object to be placed at any position and scaled to any size. 

• Field 20 is a pair of real numbers (on two lines) specifying the real and imaginary 
parts of — for the T.M case or for the TE case. 

• Field 21 is a pair of real numbers (on two lines) specifying the real and imaginary 
parts of e, for the T.M case or for the TE case. 

• Field 22 is a character flag that if set to "P" or "p" enables a plane wave generator. 
The plane wave is propagating down the y axis, and generates an condition 
on the offset boundary contour. If the flag is set to "C" or "c", a cylindrical mode 
generator is enabled. This generates an £„ cos ncf) condition on the offset boundary 
contour. If a "P" , "p", "C" or "c" is not detected, then manually input boundary 
conditions must follow, and fields 23 and 24 are not used. 

• Field 23 is a real number specifying the wave amplitude. 

• Field 24a is a real number specifying the wave frequency (in Hz). 

• Field 24b (only for cylindrical case) is an integer specifying the mode number, n. 

• Field 25 is the object perimeter data, in either polar or rectangular form. 
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An example of an input data file for a homogeneous circular cylinder is provided in 
Appendix A. The majority of the input data is echoed to the computer display and a 
system "pause" is initiated to allow for user inspection. The program may be aborted 
or continued at this time. The initial object dimensionalization has already occurred, 
and the number of unknowns and the maximum unknown width is displayed. These 
factors give an excellent indication of the expected run time for the FEBV routines. For 
example, a problem with 512 unknowns and a maximum unknown width of 31 took 806 
seconds to execute while a run with 8 unknowns and a maximum unknown width of 3 
took only 10 seconds to execute. These times are for a Intel 80386 based personal 
computer with an Intel 80287 co-processor chip calculating the fields for a circular cyl- 
inder. The FEBV routines are the next code block to execute after the pause is cleared. 

C. MESH GENERATION PROGRAM 

The mesh generation program consists of seven subroutines. These subroutines are 
an integral part of the finite element program and, therefore, were not separated. These 
routines are discussed below. 

1. lO (Input/Output) 

This subroutine reads the information contained in the 1NPUT.DAT file dis- 
cussed earlier. A two dimensional object can be described in any number of ways, 
however, for simplicity, the polar and rectangular coordinate systems are used. In either 
case, the initial assumption is that all data points are referenced to a local origin. This 
local origin can be offset by any desired amount using field 18. This offset is independ- 
ent of the entered data points or any size scaling provided by field 19. A plot label file, 
named TEXT. LB L, is created and the initial object perimeter (coded for a display in 
blue) is written to the output file, PLT.DAT. These data points describe the perimeter 
of the object. This will later prove helpful in determining the conformity of the gener- 
ated mesh to the input perimeter. All subsequent screen writes are coded for a display 
in green. Two example objects are shown in Figure 9 on page 27. These objects will 
be used throughout this chapter. The circular cylinder was generated by a separate 
computer program. The "horseshoe" shaped object was manually input. Graph paper 
was used to determine the x and y coordinates of the 28 unequally spaced perimeter 
points. 
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OBJECT PERIMETERS 

Figure 9. Typical Objects 

2. Rotate 

Tills subroutine reorders the input data points to allow for any desired bisection 
segment start and/or stop node. This subroutine is only used if manual selection of the 
bisection segment start and/or stop nodes is requested. 

3. Bound (Boundary) 

This subroutine sub-divides the object perimeter based on the mesh resolution 
specified in Field 7. The mesh resolution is the approximate length, specified in wave- 
lengths, that the user desires the perimeter to be divided into. The division of the per- 
imeter is based on linear interpolation between input data points. A new perimeter node 
is placed at each of the sub-division points. Additional input data points are necessary 
in areas of rapid change to allow for a correct object perimeter representation. 1 he ob- 
ject bisection extends from the "start node" to the "stop node". These nodes are user 
specified in Fields 16 and 17 and, in general, divide the object in half I he bisection 
should be arranged such that the object width, perpendicular to the bisection segment, 
is minimized. Thus, a long slender object should be oriented for a major axis bisection. 
Whether the major axis is oriented vertically or horizontally is of no importance. This 
is demonstrated in Figure 10 on page 28 (right side). 



BISECTION SEGMENT 



ROW 



•—OBJECT PERIMETER 

Figure 10. Bisection Segment and Row Arrangement 

For more complex objects, the bisection is not a straight line, but rather a series 
of line segments that approximate a curve. The number of perimeter nodes on the left 
side of the bisection segment must equal the number of perimeter nodes on the right side 
of the bisection segment. The bisection segment also contains this same number of 
nodes. Thus, the program must adjust the requested mesh resolution to ensure proper 
nodal spacing, 'fhis allows for a piecewise continuous segment to cross the object. 

4. Normal 

This subroutine calculates the unit normals to the object perimeter at each 
newly established perimeter node. The normal is a perpendicular constructed to the 
chord connecting the two nodes adjacent to the node for which the calculation is oc- 
curring. This perpendicular originates at the current node and is of unit length. See 
Figure 1 1 on page 29. 
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5. Nodset (Node Set) 

This subroutine uses a two sweep technique to compute the number of nodes 
on each nodal row. The rows are labeled 1 = 1, 2, 3, 1 maximum (IMX) consec- 

utively from the top of the object to the bottom. Segments normal to the perimeter 
surface, calculated in NORMAL, connect the perimeter nodes to the offset boundary 
contour. Two such completed normal segments, one on each end, complete a nodal row. 
When all rows are complete, a second contour has been established that approximates 
the object's perimeter. This contour is displaced from the perimeter by the distance 
specified in Field 8. This contour provides the boundary conditions for the FEBV rou- 
tines. The optimal selection of this offset distance will be a topic of Chapter IV. 

With the object divided into rows, each row can be further divided into equally 
spaced nodes. Based on the requested resolution and whether the objects dimensions 
are expanding or contracting, the nodal spacing is adjusted to keep the elements 
approximately the same size. Two adjacent nodal rows form an elemental row. These 
rows are given the same label, I = 1, 2, 3, ..., IMX as the upper nodal row. The nodes 
on an elemental row are numbered consecutively, starting with the left-most node. 1 his 
node is always an offset boundary contour node. When the upper nodal row is 
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numbered, the process is continued for the lower nodal row. An example of this element 
row numbering scheme is shown in Figure 12 on page 30. 




A mesh orientation attribute is set for the left and right portion of each element row. 
This attribute determines if the object has a clockwise or counter-clockwise orientation. 
See Figure 13 on page 31. Switching between the four available mesh orientations al- 
lows for a mesh that more accurately conforms to the input object perimeter without 
resulting in a disproportionate mesh structure. 
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LEFT SIDE 



RIGHT SIDE 




Figure 13, Mesh Orientation Attributes 



1 here is an additional requirement that the number of nodes on a given left or right half 
row must be only one more, equal to, or one less than the number of nodes on the ad- 
jacent left or right half row. Thus, a two sweep process is utilized. This process ensures 
that this requirement is met and that the first and last row have only two nodes. The 
second and second from last row (if present) must have three nodes. It is possible, as 
shown in Figure 14 on page 32, to generate a mesh that has only three rows. This object 
was input as a circular cylinder. It is evident that, due to a small number of elements, 
the generated mesh more accurately resembles a square. This will be a problem for cal- 
culating the scattered fields, however, the internal fields can be accurately approximated. 
In this special case, the second and second from last rows are the same. It will be shown, 
in Chapter V, that since this mesh does not closely approximate the input objects per- 
imeter the resulting Green's function contour integrals and subsequent far field calcu- 
lations have reduced accuracy. 

The nodes of the nodal rovvs form the vertices of triangular elements. An ele- 
ment, as defined in Chapter II, has three vertices, each assigned a unique number (1, 2 
or 3). The ordering of these vertices depends on the mesh orientation attribute. Each 
element also has an associated relative dielectric constant (£,) and relative permeability 

(mJ- 
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ROW 1 




Figure 14. Three Row Cylinder 



ROW 3 



6. Sorter 

This subroutine generates a complete mesh row and all the element/node inter- 
connection relationships. With the nodal structure in place, individual nodes can be 
assigned to individual elements within the global mesh structure. Each element in a 
given row is assigned a unique local element number starting with the left most element 
(relative to the bisection segment). See Figure 15. 
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It is vital that a method of determining which nodes form the vertices of a given element 
and which elements have a vertex attached to a given local node. The nodes that are 
not part of the offset contour have unknown field values. A single node can be con- 
nected to as many as six elements, only four of which can be in the current row. 1 his 
relationship is shown in Figure 16(a). This hexagonal arrangement of six elements is 
very common within the mesh. Along the bisection segment or where the mesh orien- 
tation attribute changes from one row to another, this pattern is disrupted. An example 
of an extreme case (the center of a circular cylinder) is shown in Figure 16(b). After all 
elements and nodes in an elemental row are assigned, an ordered sweep is conducted of 
that elemental row to determine which nodes are connected to which elements, which 
elements are connected to which nodes, and how many elements a single node is con- 
nected too. A complete row’ has now been generated. See Figure 17 on page 34. 





THIS NODE THIS NODE 

a. b. 

Figure 16. Two Possible Element Intersections, (a) extreme, (b) normal 

The information necessary to generate the rows and elements will be used again, in 
Chapter IV, to solve the (FEBV) problem. 

7. Finder 

This subroutine determines the x, y coordinates of each node. This data is also 
necessary to solve the FFBV program of Chapter IV and provides the plotting 
coordinates for the PLT.DAT file. This process is repeated for all rows. Thus, the entire 
object is generated as the compilation of elemental row’s made of triangles. See 
Figure 18 on page 34. 
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Figure 18. Global Mesh Structure 

A file is now available for display using a commercially available program called 
"CURVE-DIGi riZER". Any program that can accept x, y coordinate data will 
accomplish the same display process. Minor changes to the MESM program provided 
in Appendix B may be necessary since several "CURVE-DIGITIZER" plotting codes 
are embedded in the file generation code. 
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D. OPTIMIZATION OF THE MESH 

For most objects, it is recommended that the MESH program be executed in the 
mesh generation only mode (Field 6 set to "M" or "m") prior to the execution of the total 
finite element program. This will ensure that the desired mesh structure is obtained prior 
to solving for the unknown field values. The MESH (generation only) program requires 
only a few seconds for even the most dense mesh structures. It is readily seen that as 
the mesh resolution increases, the number of rows increases linearly while the number 
of unknowns increases geometrically. Therefore, the mesh calculation times may not 
change appreciably when the mesh is made more dense, however, the calculation time 
for the finite element program will rise geometrically. For this and other reasons, the 
mesh density should be kept low. This will lead to a smaller number of unknowns and 
result in faster program execution times. Figure 19 on page 36 plots this relationship 
for a circular cylinder. The solution for these unknown field values will be a topic of 
Chapter IV. 

Six different mesh generation methods are available to create mesh structures. Some 
of these methods were evolutionarx’ in nature and provide limited practical benefit. 
Methods 1 and 6 are by far the most useful. Each method is discussed below. 

Method 1 constructs the bisection segment by connecting the first input data point 
to the midpoint of the perimeter. This segment is divided into equal length segments 
each separated by a node. This method is useful for very simple objects such as the 
circular cylinder shown in Figure 20 on page 37. This circular cylinder will be used in 
the explanations of all mesh generation methods. These illustrations are not intended 
to optimize the mesh structure but rather allow for easy comparison of the 6 basic 
methods. 

.VIethod 2 constructs the bisection segment as described in method 1. The bisection 
segment nodes are. however, ordered differently. This segment is divided by connecting 
line segments from corresponding nodes on the left and right side perimeter. Where 
these line segments intersect the bisection segment, a node is placed. This leads to un- 
equally spaced bisection segments. This can be seen in Figure 21 on page 37. 
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Figure 19. Uiikiionns and Number of Rous Versus Mesh Resolution 



Method 3 modifies method 1 by allowing the user to specify, using Field 17, the stop 
node for the bisection segment. I his method can be useful to rapidly adjust around 
slightly irregular objects. This can be seen in Figure 22 on page 38. 

Method 4 combines methods 2 and 3 by using the connected line segment technique 
of method 2 to determine the node positions on the bisection segment, the stop node of 
which is specified as in method 3. This can be seen in Figure 23 on page 39. 

Method 5 improves upon method 4 by repositioning the unequally spaced bisection 
nodes. Linearly interpolated positions for the nodes leads to equal spacing. This 
method reduces the node "bunching" that frequently occurs with methods 3 and 4. This 
can be seen in Figure 24 on page 40. 
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riguie 20. Method 1 Mesh Structure Example 
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Figure 21. Method 2 Mesh Structure Example 
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Figure 22. Method 3 Mesh Structure Example 

Method 6 is the final improvement in which method 5 was modified to allow for a 
user- specified start node. This provides the user with the ability to start and stop the 
bisection segment at any input node without having to rearrange the input data. Since 
method 6 contains all of the capabilities of the other five methods, it is almost exclu- 
sively used for mesh generation. This can be seen in Figure 25 on page dO. There are 
situations that could be best served by one of the other methods. For example, a cir- 
cular cylinder and other simple symmetric objects can be represented using method 1. 
Method 2 would work well with a square or diamond shape. These objects would be 
bisected by a diagonal. Method 3 would be most suited for a synunetric object with a 
planar material interface. The more dense mesh would be used for the higher 
permittivity material. Method 4 would work well with a rhombus or other slightly 
asymmetric object. Method 5 is almost as versatile as method 6, and would work well 
in any situation where the start node is fixed. A great deal of planning is not needed in 
designing most mesh structures since they are calculated and displayed in a matter of 
seconds. 



38 




Figure 23. Method 4 Mesh Structure Example 

Iterative selection of different mesh generation parameters has proven to be the best 
technique. A summary of all mesh generation capabilities is provided in Table I. 

Appendix C contains a program called READ. FOR. This program takes the output 
data file from the "Curve-Digitizer" CAD program, called F1NALDWG.DAT, and after 
receiving the answers to several prompted questions, creates a new 1NPUT.DAT file. 
Typically, the CAD program is used to generate the perimeter of the object. This may 
be as a series of points, line segments or a combination of the two. The answers to the 
prompted questions provide the additional information needed to fill the remaining data 
fields. This is intended to be a first step towards allowing a user to specify or design an 
object and then be able to calculate the scattered fields from this object. Although it is 
far from efficient, it does serve a definite purpose. Further refinement of this program 
will allow for an easier user interface. This should enable technically trained personnel, 
without a detailed understanding of these programs, to benefit from the Field Feedback 
Formulation. 
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Figure 24. Method 5 Mesh Structure Example 




Figure 25. Method 6 Mesh Structure Example 
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Table 1. SUMMARY OF MESH GENERATION CAPABILITIES 



Method 

Number 


Straight 

Line 

Bisection 


Straight 
Line From 
Left to 
Right Side 


User Se- 
lected Stop 
Node 


Equally 

Spaced 

Bisection 

Nodes 


User Se- 
lected Start 
Node 


1 


X 






X 




2 




X 








3 


X 




X 






4 




X 


X 






5 






X 


X 




6 






X 


X 


X 
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IV. FINITE ELEMENT BOUNDARY VALUE PROGRAM 



A. INTRODUCTION 

The Finite Element Boundan.' Value (FEBV) program is the feed forward (U) oper- 
ator in the Field Feedback Formulation (F^) as shown in Figure 2 on page 4. This 
program solves the Helmholtz equation, as discussed in Chapter II, for the Dirichlet 
boundarv’ condition specified in the input data file, as discussed in Chapter III. These 
boundary conditions are imposed on the offset boundarv’ contour. The purpose of this 
program is to find the unknown field values inside the offset boundarv' contour within 
the input object. Since the ultimate goal is to obtain the scattered far fields from the 
object, the unknown field values on the perimeter are of primarv' interest. It will also 
be necessary to approximate the normal derivative of the field at the object perimeter. 
As derived in Chapter II. the goal of this program is to solve, 

+ [C,]T,+, = - [P,]T, 

The A matrix represents the effect that the I — 1 row has on the I''’ row values. Similarly, 
the B matrix represents the effect that the I-th row has on itself, while the C matrix re- 
presents the effect that the 1 + 1 row has on the I-th row. All other rows do not effect 
the I-th row. This is due to the pyramidal basis functions, discussed in Chapter II, that 
all have zero value when a node is not directly connected to another node. The P matrix 
represents the combined effects of the boundary nodes on the 1-th row. These values 
are transposed to the right side of the equality to form the system forcing function. It 
is worth remembering that for the I = 1 row, the A matrix = 0 and that for the 1 = 
I. MX row, the C matrix = 0. Thus, using the row by row stepping process, first dis- 
cussed in Chapter III of the mesh generation process, the A, B, C and P matrices can 
be filled. There is an A, B, C and P matrix for each row, and it is necessarv- to calculate 
the functionals for two elemental rows to fill one set of matrices. It is, therefore, obvious 
that the data is used twice, once for the current row and again for the ne.xt row. Spe- 
cifically, the functionals necessarv' to complete the B matrix and totally fill the C matrix 
is used again to totally fill the next row's A matrix and partially fill the B matrix. Thus, 
for a row I - 1, (such that I - 1 1), all of the A matrix and part of the B and P matrices 

are filled. When the row is stepped to I, the B and P matrices are completed and all of 
the C matrix is filled. The A, B, C and P matrices represent tri-block matrices within a 
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much larger global matrix. This row of tri-block matrices in the global matrix can be 
partially solved using the forward portion of the Ricatti transform. The Ricatti trans- 
form is a numerical technique that optimizes the solution for banded (tri-block) systems 
of linear equations. The equations are, 

Kv=-{B, + A,R,r'C, 



and 



where is the i''’ -I- 1 R matrix and the S,_, is the /'* + 1 S vector. Note that the next 
rows R matrix and S vector depends on the previous rows R matrix and S vector. 
During the forward step (MARCH subroutine) an R matrix and S vector is generated 
for each row. When all rows have been calculated, a back sweep computes the unknown 
field values, The equation is 



1 — T ^i‘ 

This back sweep must read the RS data from the disk backwards. This is a storage in- 
tensive process for all but the most trivial problem. The field values is first found 
by remembering that the last row ( 1 = IMX ). C,^,y = 0. With C;v^;^ = 0 , 
and = Ri\tx-\ • 'vhich is the last calculated S vector. The recursion continues until 
all of the \l/i 's have been calculated. 

B. FINITE ELEMENT BOUNDARY VALUE PROGRAM 

The eight subroutines comprising the FEBV program are discussed below. 

1. Zero 

This subroutine fills all the array positions of the A, B, C and P matrices with 
0 + jO. This zeroing is necessary after all calculations concerning a given row are 
completed. 

2. Varint (Variational Integration) 

This subroutine calculates the complex functionals for a given input element. 
The functionals reflect the effect that each node has on the three nodes associated with 
an element. The subroutine must be provided with the X, Y coordinates of the element 
nodes, and the material parameters e, and These functionals are returned in a 3 x 3 
complex matrix. The area of the input element is also provided as a by-product of this 
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numerical integration. The areas of each element are sorted and used to find the largest 
and smallest elements. Once found, the largest and smallest areas are combined to form 
the area ratio, which is defined as, 



area ratio = 



maximum area 
minimum area 



This ratio should be kept as small as possible. For a circular cylinder, values slightly less 
than 2.0 are possible. For more complicated objects, the ratio can get considerably 
larger. A ratio > 2.5 causes a screen message of "YOU SHOULD CONSIDER 
ABORTING THIS RUN .AND LOOKING AT THE MESH IN CURVE DIGITIZER. 
A BETTER .METHOD MAY BE AVAILABLE". It may or may not be possible to 
obtain an area ratio < 2.5. The intention of the area ratio is to insure that a uniform 
mesh is constructed prior to attempting a problem solution. Smaller area ratios are in- 
dicative of mesh structures that do not have grossly different element sizes. This will 
lead to more accurate finite element solution. 

3. Fill 

This subroutine calculates and stores all of the functionals for an entire ele- 
mental row. This is accomplished by repeated VARINT calls for each element that 
comprises an elemental row. 

4. BNDC (BoundaiT Condition) 

This subroutine calculates and stores the boundary conditions desired for the 
offset boundary contour. These conditions can be plane wave, cylindrical modes or in- 
put manually. For plane wave conditions, the percent error of the FEBV program will 
also be returned. This calculation only has meaning for the uniform case (Field 5 set to 
"U" or "u"). .Memory limitations do not allow this feature for the cylindrical mode 
boundary conditions for modes other than zero and one. A separate routine could be 
made available for offline percent error calculations. No such feature is provided if the 
boundary conditions are manually input. Any desired boundary condition may be gen- 
erated by modifying or appending the proper code to this subroutine. 

5. Loader 

This subroutine loads the A. B, C and P matrices for a given row. For all but 
the 1 = 1 and I. MX rows, two calls of FILL/ VARINT for two elemental rows of data 
are required to fill the A, B. C and P matrices. This is an additive process that starts 
after the ZERO subroutine initializes the matrices. Successive functionals are added to 
the existing data in the proper array positions. 
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6. March 



This subroutine performs the forward portion of the Riccati transform. The 
output data, called RS.DAT (R matri.x and S vector), are stored on disk. 

7. CSMINV (Complex Square Matrix Inversion) 

This subroutine accomplishes the complex matrix inversion required by the 
forward Riccati transform. A maximum dimension of 50 x 50 was established to limit 
memory utilization. This allows for a maximum unknown width of 50. 

8. Sweep 

All of the above subroutines are called at least once for each row. The sweep 
routine is called only after all of the row calculations are completed, and then only once. 
This subroutine conducts the Riccati back sweep by reading the RS data generated by 
the MARCH subroutine. This data is read olf the disk backwards, using the 
FORTRAN "backspace" command. The returned field values are actually individual 
contributions due to a unit valued basis function being individually applied to each of 
the boundary nodes. In so doing, the problem need only be solved once. After the data 
is stored, in matrix form as the U.DAT file, the boundary value problem may be solved 
for any incident field by multiplying the U matrix by the new incident fields. Thus. 



^perimeier Cff] * ^ incidenr 



A summar\‘ of the input, output and error data is provided in the OUTPUT, DAT file. 

9. Save 

This subroutine was necessar\‘ to allow for reasonably sized problems. Since 
all of the Field Feedback Formulation code could not fit into 640 kilobytes of memor>', 
the program was divided into two parts. All of the data necessary to perform the pro- 
grams is saved in the F3.DAT file. This data is then read by the field feedback program. 
This technique, though very inefficient, circumvents the 640 kilobyte memor}’ limitation 
imposed by the IB.M Disk Operating System (DOS). Efforts to convert this code to run 
under a compiler that does not have a 640 kilobyte memory^ limitation, such as Micro- 
way NDP FORTR.AN compiler, will be pursued at a later date. 



C. VARINT VALIDATION 

The first step in the program validation required an understanding of the error 
convergence of the variational elements as a function of element size, d and material 
properties. £, and a, . The test program provided in Appendix D varies the element size, 
in wavelengths, for a test mesh structure. This test structure shown in Figure 26 on 
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page 46 has only one unknown at the center. The mesh is made of four adjacent ele- 
ments. These elements are inside a uniform material having properties, e, and /j,. The 
other four nodes are established as boundary nodes. 

PLANEWAVE 



Figuie 26. 




A plane wave, of user specified frequency and amplitude, is used to determine the 
boundary conditions on the four boundary nodes. This plane wave is propagating down 
the vertically oriented axis. The unknown nodal value is calculated and compared to the 
actual plane wave value (1 -t- jO). A percent error is calculated as the dimension of the 
elements, d are reduced. As expected, the solution became more accurate as the ele- 
ments become smaller. A plot of the results of this test program is provided in 
Tigure 27 on page 47. This data was generated with c, = 1 -f JO, and n, = I + JO. 
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Figure 27. Solution Error for a Test Mesh Structure 

The region of interest is where the error approaches zero. An expanded plot of this re- 
gion is provided in Figure 28 on page 48. The accuracy of the solutions were desired 
within 1 percent of the exact value, 'fhis requires an element that is < -jy , in the ma- 
terial. To ensure the desired error was achieved, elements were usually scaled to be 
< , in the material. The phase of the plane wave was also varied to determine the 

effect on convergence. No significant effect was noted. 
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Figure 28. 



Solution Error for a Test Mesli Structure (Expanded) 



D. FINITE ELEMENT BOUNDARY VALUE PROGRAM VALIDATION 

1 he next validation step required an actual object to calculate the fields in and on. 
The simplest object was actually not an object at all, but rather an imaginary circular 
cylinder in free space. A plane wave was propagated through this free space. Since no 
material interface existed, the exact solution would be known for all positions. The 
plane wave established the boundary conditions on the offset contour. Trial runs were 
conducted varying the mesh resolution and boundary offset contour distance. With the 
error defined as. 



error = 



^(calculated value — actual value)^ 



I(actual value)^ 



two errors were defined. The perimeter error only considers the perimeter nodes in the 
error calculation. The bisection segment error only considers the bisection segment 
nodes, with the exception of the two end nodes, in the error calculation. The end nodes 



48 



of the bisection segment are part of the perimeter and are, therefore, not considered. 
As expected, tlie perimeter error could be rapidly reduced by decreasing the olfset dis- 
tance. 1 his is due to the fact that the node or nodes closest to an unknown node dom- 
inate the field contributions at this node. Again, a goal of < 1% error was desired. 1 he 
reduced offset distance had a very small effect on the bisection segment error. The only 
way to significantly reduce this error was to increase the mesh resolution. An increased 
mesh resolution reduced both errors. Figure 29 show's the perimeter error as a function 
of the number of bisection segment nodes. 




Figure 29. Perimeter Error 

The three curves are for offset distances of 0.05, 0.025 and 0.01 X . Figure 30 on page 
50 show's the bisection segment error as a function of the number of bisection segment 
nodes. The two curves are for offset distances of 0.05 and 0.025 . For offset distances 

less than 0.025 2o , the curves are almost identical to the 0.025 curve. These curves 
are omitted for clarity. 
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Number of BltooHon Nodot 

Figure 30. Bisection Segment Error 



The calculated and exact field values for a 0.5 diameter circular cylinder, with a 
0.05.^,, mesh resolution and ofi'set distance and c, = 1 + jO is shown in Figure 31 on page 
51. The exact fields values for the real and imaginary portion of the plane wave are 
shown as squares and diamonds, respectively. 1 he solid curves are the calculated field 
values. The perimeter and bisection errors were 0.74 and 1.69 percent, respectively. 
Figure 32 on page 52 shows the efiects of not properly adjusting the mesh resolution 
and offset distance when the material is changed. In this case, the permittivity was 
changed to e, = 4 + jO. 1 his caused the perimeter and bisection segment errors to in- 
crease to 18.2 and 41.1 percent, respectively. The mesh resolution and olfset distance 
were reduced to 0.02 U.q and the permittivity was changed to £, = 4 — J4. The results are 
shown in Figure 33 on page 53 . This proper selection of mesh resolution and olfset 
distance reduced the perimeter error and bisection segment error to 0.44 and 0.56 per- 
cent, respectively. Note that in the lossy material case the two errors are very close in 
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magnitude. This is due to the way that the error is defined. This definition, in a lossy 
material, overemphasizes the objects leading edge. 




Figure 31. 0.5 ). cylinder, c, = 1 + jO 



This portion of the object (as well as the trailing edge) is the most accurate for the 
bisection segment calculations. This is because of the relative closeness of these nodes 
to the perimeter. For the lossless material, the bisection segment error is typically two 
to three times the perimeter error, for the same number of bisection nodes. The con- 
clusion that the oflset distance should be made as small as possible in order to minimize 
the perimeter error is correct, but will prove to be counterproductive in the long run. 
There is a competing effect that will require this contour offset distance to be as large 
as possible. This effect, associated with the accuracy of the Green's function contour 
integral, will be discussed in Chapter V. 
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Figure 32. 



0.5 A cylinder, £, = 4 + jO 
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Figure 33. 0.5 ). cylinder, £, = 4 — j4 



E. INHOMOGENEITY 

Until this point, all testing was done in circular homogeneous dielectric materials. 
It was, at a minimum, desirable to place the object in a vacuum to calculate the scattered 
fields. This requires the first two and last two elements of each row to have 
E, = 1 + yO and = 1 + yO. These elements represent the space between the objects per- 
imeter and the ofTset boundary contour. Since the plane wave solution is no longer valid 
for this geometry, a new problem with a known solution was needed. 

Given a homogeneous dielectric circular cylinder of radius a and permittivity e, , as 
shown in Figure 34 on page 54, it can be shown that [Ref 8 ], 

4>) = A„J„{k,R) cos n(f), 
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where, 



and 



ii„{R, <f>) is the field value inside the dielectric material. 



_ -j2D„ 

^ Ji,{R,)[j„{k,R,)U'^hR,) - krJ' „{WH'^\r,)] - 
' H^^\Rt)lUKRaVn{K) - k^\{KR,)JM-] 

J„ = Bessel Function of the kind of order n 



k, = J7,. 
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An additional formulation is required for the fields outside of the cylinder. A cylindrical 
wave (mode, n = 1) was applied to a homogeneous dielectric cylinder. The exact and 
finite element field values were calculated and compared. Figure 35 on page 55 , shows 
a typical result. Mesh resolution, olfset distance and material permittivity were varied 
to determine convergence. The results were similar to the dielectric homogeneous plane 
wave case discussed earlier. 




P«r1m«t«r Nod« Number 

Figure 35. Typical Cylindrical Mode Boundary Value Problem Result 

F. FINITE ELEMENT CONCLUSIONS 

High accuracy solutions can be obtained by maintaining a mesh resolution of 
<-^ in the material. The offset distance should be kept at approximately the same 
magnitude as the mesh resolution, but may be as large as . Increased accuracy re- 
quires an increased mesh resolution. This increase in accuracy is at the expense of in- 
creased processor time. 
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V. FIELD FEEDBACK PROGRAM 



A. INTRODUCTION 

The Field Feedback Program is the feedback (T) operator in the This program 
utilizes the output of the finite element program and the input incident fields to evaluate 
a "near field" Green's function integral. This integral is called "near field" in that the 
integral will be performed within of the object perimeter. As seen in Figure 36, 
three surface contours are defined. The boundary contour is where the incident field 
boundary condition is applied. 




Figure 36. Contour Arrangement 



The perimeter contour is where the finite element program solves for the field values on 
the object perimeter. The geometric contour is midway between the boundary and the 
perimeter and is the contour over which the Green's function contour integral is per- 

C'A . • -rr- 

formed. This is necessary to allow the to be approximated by a finite diflerence 
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technique. The i/^^c is the average of the field values on the boundary and perimeter. 
Thus, 



O^j/ boundary ^ perimelei^ 

on dn 



and 



•Acc — 



perimewr ^boundary) 



The Field Feedback program is provided as Appendix E. 

B. FIELD FEEDBACK PROGRAM 

1. Input 

This subroutine reads the finite element data stored in the F3.DAT file by the 
S.AVE subroutine. All necessary nodal X. Y coordinates are calculated. 

2. TMAT (T Matrix) 

This subroutine loads a complex matrix called T.MAT. Each element of the T 
matrix is the result of a Green's function contour integral. This integral is conducted 
on the GC contour. The m-th column of the T matrix is evaluated with a single unit 
valued basis function on the m-th boundary node. The equation, 



T = 

^ n, m 






^GC 



dGC, 



on 



on 



may now be numerically evaluated at each of the n boundary' nodes. This process is 
repeated for each row until the entire matrix is filled. 

The numeric integration uses a rectangular mid-point approximation technique. 
For this linearized problem, this is equivalent to a trapezoidal integral technique. The 
integral path between any given two nodes is subdivided based on the distance to the 
first point. This subdivision maybe controlled by three input variable fields. A multi- 
plicative scale factor and offset terms are available. To date, the utilization of the scale 
factor (set > 1.0) and the offset terms (set > 0 ) have not proved necessar}'. This may 
be necessary for more complicated geometric structures. 
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3. CNSOLV 



This subroutine solves the equation, 



matrix^ * ^ boundary^ 



where C„ is the total scattered fields on the boundar}', I is the identity matrix and 
^boundary is thc initially specified incident fields. 

4. FFLD (Far Fields) 

This subroutine calculates the scattered far field or bistatic radar cross section 
per unit length of the object. Any desired integer angular resolution for the calculation 
may be specified. The final results are stored in the FFPAT.DAT file. 



C. FIELD FEEDBACK VALIDATION 

Given a plane wave boundary conditions imposed on an imaginary- cylinder in free 
space, a Green's function contour integral may be performed around the cylinder. This 
cylinder is imaginary in the sense that it does not actually exist. The cylinder is actually 
an artificial boundan.' that does not form a material interface. For points inside this 
cylinder, the resulting integral should equal the negative of the actual plane wave value 
at that point in freespace. For points outside the cylinder, the resulting integral should 
equal zero. Several test conditions were investigated to ensure that the numerical inte- 
gration of the Green's function did arrive at the expected values. Maximum absolute 
errors for these tests were less than 7 x 10"'’ . These test cases started with exact values 

Cl// 

of the field. \p and the normal derivative, — . After initial validation, finite element 

Cl/' 

calculated xh and — values were used. Errors increased bv approximated a factor of 

Cn . rr . 

10. Numerous competing effects were observed with two dominant effects becoming 
evident. 



1. Small Object Phenomenon 

For very small objects, where only a few elements are needed to meet the 
mesh resolution requirement, the normal derivative accuracy is crucial. Since the normal 
derivative is the calculated normal to the piecewise linear approximation of the objects 
perimeter, this fit is usually not adequate. To prevent this problem, additional 
perimeter; boundary nodes are needed. This is easily accomplished by increasing the 
mesh density, which is controlled by the mesh resolution (input field 7). For these ver>‘ 
small objects, the RCS is almost uniform for the T.M case and co-sinusoidal for the TE 
case. The utility of this calculation must, therefore, be questioned. 
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2. Offset Distance Phenomenon 

The ofTset distance must not be made too small or the accuracy of the Green's 
function integral will suffer. For a circular cylinder, the optimal offset was between 0.03 
}. and 0.035 /. Note that these distances are always < . 
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VI. VALIDATION 



A. INTRODUCTION 

The difTiculty in the total Field Feedback Formulation (F^) program validation is in 
finding problems that have well established or accepted solutions. The total program 
is the combination of the mesh generation 'finite element program and the field feedback 
program. If possible, a problem that has a closed form analytical solution is desired. 



B. HOMOGENEOUS CIRCULAR CYLINDRICAL SCATTERING 

This is the first test case for the program. A penetrable circular cylinder has a 
closed form analytic solution, so the final results could be verified against exact values. 
It can be shown that for the TM case (E - wave) the reflection coefficient is. 




UWJ'niR,) - J'JWUR,} 

MWIi'FIR.) - J'n(KR,)ii°\R,) 



and that for the TE case (H - wave) the refection coefficient is, 



F 



TM 

n 






where k, — r, = (Ref. 9 ] The scattered field for the 

TM case is. 



4(r, </>) = -£' 



incident 



Fo//n^) + 2 



oo 

z- 



r„//; 



(2) 



\R) cos n4> 



n=l 



Similarly, the scattered field for the TE case is, 



Hz(R, 4>) = -H 



incident 



d^\R) + 2 






CO 

z- 



^H^^\R) cos n4> 



n=l 
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In the far field region, as R approaches infinity, the scattering width may be defined as, 



o{4)) = 2n lim p 

p-*oo 



l/7'P 



I ^incident | 2 



and 



o{4>) = ^ |ro + 2^r„ cos«0|l 

The source code, without the Bessel function routines, for these calculations is provided 
as Appendix F. 

Three different radius dielectric cylinders were tested. Wavenumber normalized radii 
of 0.5, 1.0 and 2.0 with z, = 2.56 + jO were selected. The TE and T.M results of the 0.5, 
1.0 and 2.0 radii problems are provided as Figure 37 on page 62, Figure 38 on page 
63, and Figure 39 on page 64 . To validate the Chapter VI conclusion that the unit 
normal was the cause of the errors for ver\' small circular cylinders, the c, was increased 
from 2.56 to 25.6. This value still meets the requirement for mesh resolution not to ex- 
ceed . The TE and TM results are provided as Figure 40 on page 65 . Actual cross 
section values are indicated by the squares and the Field Feedback Formulation sol- 
utions are plotted as a single solid curve. 

C. HOMOGENEOUS IRREGULAR OBJECTS 

The results detailed in [Ref 10] were next validated with the program. This ob- 
ject, as seen in Figure 41 on page 66 , is a dielectric shell with inner radius = .25 Zq, 
outer radius = .30 and e, = 4 + jO. The comparison of the [Ref 10] data and the 
F^results are provided as Figure 42 on page 67. 
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Figure 37. Cylinder, TE and TM Case, k^a = 0.5, e, = 2.56 
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Figure 38. Cylinder, TE and TM Case, k^a = 1,0, f, = 2.56 
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Figure 39. Cylinder, TE and TM Case, k^a = 2.0, z, = 2.56 
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Figure 40. Cylinder, TE and TM Case, k^a = 0.5, e, = 25.6 
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PLANEWAVE 



V 




Figure 41. Dielectric Shell Mesh 

Tlie TE case shows excellent agreement. The TM case tends to diverge at the smaller 
values of 4>. 

D. AN INHOMOGENEOUS OBJECT 

The results detailed in (Ref. 11] were next validated with the program. This ob- 
ject, as seen in Figure 43 on page 68, is a dielectric ring with inner radius = .25 , outer 

radius = .30 and e, = 4 + jO. The exact solution to this problem is available. [Ref 

11 ] 
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Figure 42, Dielectric Shell, TE and TM Case, c, = 4 + jO 
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Figure 43. Dielectric Ring 

The mesh generation program was not modified to conform to both the outer and inner 
material boundaries. The original mesh generation program design concept was to 
closely fit the mesh to the objects outer perimeter and then allow for material inhomo- 
geneities to be accounted for by difi'erent material parameters being assigned to individ- 
ual elements. A section of the generated mesh is shown in Figure 44 on page 69. Note 
the additional curved line showing the .25 2o inner radius. This curve does not follow 
the established element boundaries. The efiective ring is shown in Figure 45 on page 
69. This resulted in the inner radius having an irregular pattern as elements near the 
material interface varied in material composition. This is similar to the granular noise 
problem characteristic of delta modulation communication channels. The TE and TM 
results are provided as Figure 46 on page 70. 
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Figure 45. Effective Geometry for the Dielectric Mesh 
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Figure 46. Dielectric Ring, TE and TM Case, e, = 4 + jO 
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VII. CONCLUSIONS 



A. RESULTS 

The Field Feedback Formulation proved to be an excellent tool for calculating the 
internal and scattered fields from the tested inhomogeneous asymmetric objects. The 
keys to satisfactor}’ results are, 

• Keep the maximum element dimension < , 

• Maintain a mesh with a uniform distribution, and 

• Maintain the offset distance near the mesh resolution in magnitude, but no greater 
than . 

The techniques used proved to be capable of adapting to a wide variety of situations. 
These adaptations did not require program modifications or reprogramming. The nu- 
merous built-in features, as discussed in the input field description of Chapter 111, al- 
lowed for this robustness. 

B. RECOMMENDATIONS AND EXTENSIONS 

With the basic program validated, future testing and development should, 

• Emphasize large object validation against accepted results. 

• Stress inhomogeneity and irregularity in all testing. 

• Optimize the T matrix element calculation by improving the Green's function 
contour integral. 

• Evaluate the usefulness of the maximum distance feature (Field 12). Beyond this 
maximum distance no contribution is made to the Green's function integral. It is 
not expected that this will be possible for objects less than a few wavelengths in 
maximum dimension. 

• Modify the mesh generation program to allow for conformity to multiple interfaces 
(multi-layered objects without the granular noise problem). 

• Modify all programs for Intel 80386 based Fortran compiler use, such as NDP 
FORTRAN by Microway. This will remove the 640 KB memory limitation im- 
posed by the IB.M DOS and allows for the solution to larger scattering problems. 
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APPENDIX A. INPUT DATA FILE EXAMPLE 



CIRCLE 

R 

NI 

ND 

NU 

NM 

0. 019 

0. 03 

1 . 0 
0 

0 

999. 9 
36 
72 
1 
5 

35 
4. 5 
4. 0 

0. 3 

1 . 0 
0 . 0 
1 . 0 
0 . 0 
P 

1 . 0 

3. OOE+08 



PLOT TITLE 

RECTANGULAR INPUT DATA 
NO INTERMEDIATE DATA RECORDING 
NO DISSPLA PROGRAM GENERATION 

NON UNIFORM MATERIAL (MATERIAL INTERFACE PRESENT) 

NO MESH GENERATION ONLY 

REQUESTED MESH RESOLUTION (IN WAVELENGTHS) 

CONTOUR DISTANCE (IN WAVELENGTHS) 

MULTIPLICATIVE GFI SCALE FACTOR 
DISTANCE <1.0 BIAS TERM 
DISTANCE > 1.0 BIAS TERM 

MAX DIST BEYOND WHICH THERE IS NO CONTRIBUTION TO GFI 
NUMBER OF INPUT DATA POINTS 

NUMBER OF POINTS FOR SIGMA CALCULATION (IN THE CIRCLE) 

ELEMENT GENERATION METHOD 

START NODE 

STOP NODE 

X AXIS OFFSET 

Y AXIS OFFSET 

DESIRED DIMENSION (ORIGIN TO FIRST POINT , WAVELENGTHS) 

REAL PART OF ALPHA 

IMAGINARY PART OF ALPHA 

REAL PART OF BETA 

IMAGINARY PART OF BETA 

PLANEWAVE GENERATOR ENABLED 

PLANEWAVE AMPLITUDE 

INPUT FREQUENCY OF THE PLANEWAVE, FO 

- ANGLE, X, Y COORD. 



0 


. 00000000 


. 50000000 


10 


. 08682409 


. 49240390 


20 


. 17101010 


. 46984630 


30 


. 25000000 


. 43301270 


40 


. 32139380 


. 38302220 


50 


. 38302220 


. 32139380 


60 


. 43301270 


. 25000000 


70 


. 46984630 


. 17101010 


80 


. 49240390 


. 08682407 


90 


. 50000000 


-. 00000002 


100 


. 49240390 


-. 08682411 


110 


. 46984630 


-. 17101010 


120 


. 43301270 


-. 25000000 


130 


. 38302220 


-. 32139380 


140 


. 32139380 


-. 38302220 


150 


. 25000000 


-. 43301270 


160 


. 17101000 


-. 46984630 


170 


. 08682404 


-. 49240390 


180 


-. 00000004 


-. 50000000 


190 


-. 08682413 


-. 49240390 


200 


-. 17101010 


-. 46984630 


210 


-. 25000000 


-. 43301270 


220 


-. 32139380 


-. 38302220 


230 


-. 38302220 


-. 32139380 
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240 

250 

260 

270 

280 

290 

300 

310 

320 

330 

340 

350 



43301270 

46984630 

49240390 

50000000 

49240390 

46984630 

43301270 

38302220 

32139380 

24999990 

17101000 

08682401 



25000000 
17101000 
08682403 
. 00000007 
. 08682416 
. 17101010 
. 25000010 
. 32139390 
. 38302230 
. 43301280 
. 46984630 
.49240390 
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APPENDIX B, READ PROGRAM 



c 



c 



c 



REAL A(0:2000), B(0:2000), MRES, DIST, XORIGIN, YORIGIN 
REAL C, D, E, F, DPER, DD, FREQ, EO, MAXD 

INTEGER I, J, METHOD, STARTND, STOPND, NRES, MODE, LBIAS, GBIAS 
CHARACTER*! CHAR , CHARI , CHAR2 , CHAR3 , CHAR4 , CHARS 
CHARACTER* 13 NAME 

OPEN(UNIT = 1, FILE = 'D: FINALDWG.DAT’) 

OPENCUNIT = 2, FILE = 'C: MSFORT INPUT.DAT’) 

J = 0 



WRITE(*,*) 'ENTER THE PLOT NAME OR LABEL (MAX OF 13 CHARACTERS)’ 
READ(*,1005) NAME 

WRITE(*,*) 'ENTER THE COORDINATE SYSTEM IN USE. (R OR P)' 
READ(*,1000) CHAR 

WRITE(*,*) 'DO YOU WANT INTERMEDIATE VALUES ? (I)' 

READ(*,1000) CHARI 

WRITE(*,*) 'DO YOU WANT A DISSPLA PROGRAM GENERATED ? (D)' 
READ(*,1000) CHAR2 

WRITE(*,*) 'UNIFORM SLAB (NO MATERIAL TO VACUUM INTERFACE ? (U)' 
READ(*,1000) CHAR3 

WRITE(*,*) 'MESH GENERATION ONLY ? (M)' 

READ(*,1000) CHAR4 

WRITE(*,*) 'ENTER THE MESH RESOLUTION. (0.4, ETC...)' 

READ(*,*) MRES 

WRITE(*,*) 'ENTER THE CONTOUR DISTANCE. (0.3, ETC...)' 

READ(*,*) DIST 

WRITE(*,*) 'ENTER THE GFI SCALE FACTOR. (1.0, ETC...)’ 

READ(*,*) DPER 

WRITE(*,*) 'ENTER THE GFI ( < 1 BIAS TERM (0,1, ETC...)' 

READ(*,*) LBIAS 

WRITE(*,*) 'ENTER THE GFI ( > 1 BIAS TERM (0,1, ETC...)' 

READ(*,*) GBIAS 

WRITE(*,*) 'ENTER THE GFI MAX DISTANCE (999.0, ETC...)’ 

READ(*,*) MAXD 

WRITE(*,*) 'ENTER THE NUMBER OF POINTS FOR SIGMA CALC. ( 36 ,ETC. . . ) ' 
READ(*,*) NRES 

WRITE(*,*) 'ENTER THE DRAWING METHOD. (1 - 6)' 

READ(*,*) METHOD 

WRITE(*,*) 'ENTER THE STARTING NODE NUMBER. (MUST BE > 0)' 

READ(*,*) STARTND 

WRITE(*,*) 'ENTER THE BISECTION STOPPING NODE NUMBER. ' 

READ(*,*) STOPND 

WRITE(*,*) 'DESIRED DISTANCE FROM ORIGIN TO POINT 1 (IN WLs) ' 
READ(*,*) DD 

WRITE(*,*) 'ENTER THE X AXIS ORIGIN. ' 

READ(*,*) XORIGIN 

WRITE(*,*) 'ENTER THE Y AXIS ORIGIN. ' 

READ(*,*) YORIGIN 

WRITE(*,*) 'ENTER THE REAL COMPONENT OF ALPHA. ’ 

READ(*,*) C 
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WRITE(*,*) 'ENTER THE IMAGINARY COMPONENT OF ALPHA. ’ 
READ(*,*) D 

WRITE(*,*) 'ENTER THE REAL COMPONENT OF BETA. ' 
READ(*,*) E 

WRITE(*,*) 'ENTER THE IMAGINARY COMPONENT OF BETA. ' 
READ( F 

WRITE(*,*) 'PLANE WAVE (P) OR CYLINDRICAL MODE (C) ' 
READ(*,1000) CHARS 

WRITE(*,*) 'ENTER THE WAVE FREQUENCY (IN HERTZ)' 
READ(*,*) FREQ 

WRITE(*,*) 'ENTER THE WAVE AMPLITUDE ' 

READ(*,*) EO 

IF((CHAR5.EQ. 'C').OR. (CHARS. EQ. 'c')) THEN 
WRITE(*,*) 'ENTER THE MODE NUMBER ' 

READ(*,*) MODE 

END IF 

DO 10, 1=0, 1999 

READ(1,*,ERR = 20) A(J), B(J) 

IF(A(J). GT. 1000) THEN 
GOTO S 

ELSEIF(B(J). GT. 1000) THEN 
GOTO S 

ELSE 

J = J + 1 

END IF 

S CONTINUE 

10 CONTINUE 

20 J = J - 1 

WRITE(2,100S) NAME 
WRITE(2,1000) CHAR 
WRITE(2,1000) CHARI 
WRITE(2,1000) CHAR2 
WRITE(2,1000) CHARS 
WRITE(2,1000) CHAR4 
WRITE(2,1010) MRES 
WRITE(2,1010) DIST 
WRITE(2,1010) DPER 
WRITE(2,1020) LBIAS 
WRITE(2,1020) GBIAS 
WRITE (2, 1040) MAXD 
WRITE(2,1020) J 
WRITE(2,1020) NRES 
WRITE(2,1020) METHOD 
WRITE(2,1020) STARTND 
WRITE(2,1020) STOPND 
WRITE(2,1010) XORIGIN 
WRITE(2,1010) YORIGIN 
WRITE(2,1010) DD 
WRITE(2,1010) C 
WRITE(2,1010) D 
WRITE(2,1010) E 
WRITE(2,1010) F 
WRITE(2,1000) CHARS 
WRITE(2,1010) EO 

IF((CHARS. EQ. 'C’).OR. (CHARS. EQ. 'c')) THEN 
WRITE(2,1020) MODE 
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ENDIF 

WRITE (2,*) FREQ 
DO 30, I = 1, J 

WRITE(2,1030) I, A(I), B(I) 
30 CONTINUE 

WRITE(2,*) 'END OF FILE’ 

C 

1000 FORMAT(Al) 

1005 FORMAT(A13) 

1010 FORMAT(F8. 6) 

1020 F0RMAT(I3) 

1030 F0RMAT(1X,I3,7X,F12. 8,5X,F12. 8) 

1040 FORMAT(F8. 3) 

C 

CLOSE(l) 

CLOSE(2) 

STOP 

END 

C 
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APPENDIX C. MESH GENERATION/FINITE ELEMENT PROGRAM 



FINITE ELEMENT MESH GENERATION PROGRAM 
WRITTEN BY LT. T. B. WELCH 

w/ PROGRAMMING IDEAS FROM PROF. M. A. MORGAN 
COMPLETELY FILE (INPUT.DAT) DRIVEN 
INPUT PARAMETERS: 

CHARACTER (POLAR OR RECTANGULAR INPUT DATA - FLAG) 
CHARACTER (INTERMEDIATE MATRICES WRITE - FLAG) 

CHARACTER (DISSPLA PROGRAM GENERATION - FLAG) 

CHARACTER (UNIFORM MATERIAL - FLAG) 

CHARACTER (MESH GENERATION ONLY - FLAG) 

MESH RESOLUTION 
CONTOUR OFFSET 

DESIRED SCALE FACTOR FOR GREEN'S FUNCTION INTEGRAL 
BIAS (SHIFT) FOR GFI STEP WHEN <1.0 
BIAS (SHIFT) FOR GFI STEP WHEN >1.0 
MAXIMUM DISTANCE BEYOND WHICH THERE WILL BE NO 

CONTRIBUTION TO THE GREEN'S FUNCTION INTEGRAL 
NUMBER OF POINTS 

NUMBER OF POINTS FOR CROSS SECTION, (IN THE CIRCLE) 

MESH GENERATION TECHNIQUE 
START NODE 
STOP NODE 
X ORIGIN 
Y ORIGIN 

DESIRED DISTANCE (FROM ORIGIN TO FIRST POINT, WAVELENGTH) 
ALPHA (INPUT AS A AND B THEN CONVERTED TO COMPLEX ALPHA) 
BETA (INPUT AS A AND B THEN CONVERTED TO COMPLEX BETA) 
CHARACTER (PLANEWAVE GENERATOR - FLAG) 

AMPLITUDE 

FREQUENCY (IN HERTZ) 

-- OR -- 

AMPLITUDE 

CYLINDRICAL MODE NUMBER (GENERATOR ALWAYS DOES N=l) 
FREQUENCY (IN HERTZ) 

-- OR -- 

BOUNDARY CONDITIONS 

DATA POINT PAIRS (X, Y OR RADIUS, THETA) 

END OF FILE MARKER 



COMPLEX ALPHA, BETA, BCOND( 100) ,LINE(50) ,ANS(100) 

COMPLEX A(50,50) ,B(50,50) ,C(50,50) ,P(50, 100) ,SURBC( 100) 

COMPLEX U( 100 , 100) , PSI(50 , 100) , SVEC(50 , 100) 

REAL MRES,MESH(0: 200,5) ,NDP(200,2) , OFFSET, DPER,MRESW 
REAL MINAREA , MAXAREA , AREA , EO , XORIGIN , YORIGIN , KO , MAXD 
I NTEGER NPNTS , PERND , NODE S ( 2 0 0 ) , MOR( 2 0 0 ) , B I ND , NBMX , UNK , LB I AS , GB I A S 
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c 



c 



c 



c 



c 



INTEGER LND( 0:200,3), NDL( 200,4), NCT( 200 ) , MAXEL , IMX , MODE , NRES 
INTEGER METHOD , STARTND , STOPND ,MINROW , MINED, MAXROW ,NABC( 100,3) 
CHARACTER*! CHAR, CHARI, CHAR2, CHAR3, CHAR4, CHARS 
EQUIVALENCE (PSI, P) 



COMMON/BLKl/MESH 

C0MM0N/BLK2/PERND , BIND 

C0MM0N/BLK3/N0DES 

C0MM0N/BLK4/LND,NDL,NCT 

C0MM0N/BLK5/NDP 

C0MM0N/BLK6/M0R 

COMMON/ BLK7/MINAREA , MINROW , MINED , MAXAREA , MAXROW , MAXEL , AREA 

C0MM0N/BLK8/A,B,C,P 

C0MM0N/BLK9/CHAR2, CHAR3, CHAR4 



OPEN (UNIT = 

OPEN (UNIT = 

OPEN (UNIT = 

OPEN (UNIT = 

OPEN (UNIT = 

CACCESS=’ SEQUENTIAL' ,FORM^ 
OPEN (UNIT = 8, FILE = 

CACCESS=' SEQUENTIAL' ,FORM^ 



1, FILE = 

2, FILE = 

3, FILE = 

4, FILE = 
7, FILE = 



OPEN (UNIT 
OPEN (UNIT 
OPEN (UNIT 
OPEN (UNIT 
OPEN (UNIT 
OPEN (UNIT 
OPEN (UNIT 
OPEN (UNIT 
OPEN (UNIT 



9, FILE 

10, FILE 

11, FILE 

12, FILE 

13, FILE 
19, FILE 

FILE 
FILE 



20 , 

30, 



= 40, FILE = 



MATRIX. DAT' , STATUS = 'UNKNOWN') 
PLT. DAT' , STATUS = 'UNKNOWN') 
INPUT. DAT' , STATUS = 'OLD') 

TEXT. LBL' , STATUS = 'UNKNOWN') 
RS. DAT' , STATUS = 'UNKNOWN', 
^'UNFORMATTED' ) 

PSI. DAT' , STATUS = 'UNKNOWN', 
^'FORMATTED' ) 

ABCP.DAT' , STATUS = 'UNKNOWN') 
BCOND.DAT', STATUS = 'UNKNOWN') 
FMATRIX. DAT' , STATUS=' UNKNOWN' ) 
NABC. DAT' ,STATUS=' UNKNOWN' ) 
OUTPUT. DAT' , STATUS=' UNKNOWN ' ) 
ABRMAT. DAT' , STATUS= ' UNKNOWN ' ) 
DISPLA. DAT' ,STATUS=' UNKNOWN' ) 

U. DAT' ,STATUS=' UNKNOWN' ) 

F3. DAT' ,STATUS=' UNKNOWN' ) 



MINAREA = 999999. 9 
MAXAREA = -999999. 9 



CALL 1 0( NPNTS , MRE S , METHOD , STARTND , STOPND , OFFSET , ALPHA , BETA , BCOND , 
CEO , CHARI , MODE , XORIGIN , YORIGIN , CHARS , DPER , KO , NRES , MRESW , LBIAS , GBIAS 
C,MAXD) 

CALL ROTATE ( NPNTS , METHOD, STARTND, STOPND) 

CALL B OUND( MRE S , NPNTS , METHOD, STOPND) 

CALL NORMAL 

CALL NODSET( METHOD, NABC, NBMX,UNK) 

PAUSE 'PLEASE PRESS ENTER TO CONTINUE, OR CTRL BREAK TO ABORT! ' 
CALL LOADER( BCOND , OFFSET , ALPHA , BETA , NABC , IMX , NBMX , EO , SURBC , CHARI , 
CLINE , MODE , XORIGIN , YORIGIN, SVEC , CHARS ) 

CALL SWEEP( IMX , NABC , SURBC , LINE , CHARI , U , BCOND , PS I , ANS , CHARS ) 

CALL SAVE ( BCOND , ANS , U , OFFSET , PERND , CHARS , DPER , KO , XOR I GIN , YORI GIN , 
CNRES , MRESW , LB I AS , GBIAS ,MAXD) 

CLOSE(l) 

CL0SE(2) 

CL0SE(3) 

CL0SE(4) 

CL0SE(7) 

CL0SE(8) 
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CL0SE(9) 

CLOSE(IO) 

CLOSE(ll) 

CLOSE(12) 

CLOSE(13) 

CLOSE(19) 

CLOSE(20) 

CLOSE(30) 

CLOSE(40) 

STOP 

END 

SUBROUTINE IO( NPNTS , MRES , METHOD , STARTND , STOPND , DIST , ALPHA , BETA , 
CBCOND , EO , CHARI , MODE , XORIGIN , YORIGIN , CHARS , DPER , KO , NRES , MRESW , 
CLBIAS,GBIAS,MAXD) 

THIS SUBROUTINE READS IN THE INPUT PARAMETERS 
AND SURFACE DATA POINTS. THESE POINTS CAN BE 
IN EITHER POLAR OR RECTANGULAR FORM. 

COMPLEX ALPHA, BETA, BCOND( 100) 

REAL A , B , PI , XORIGIN , YORIGIN , EO , KO , FO , SFAC , DD 

REAL MESH(0: 200,5) , MRES, THETA, RADIUS, DIST, C, LAMBDA, MRESW 

REAL DISTW,DPER,MAXD 

INTEGER 1 , 1 1 , J , K , NPNTS , RE S , METHOD , STARTND , STOPND , MODE , NRE S , LB I AS 
INTEGER GBIAS 

CHARACTER*! CHAR , CHARI , CHAR2 , CHAR3 , CHAR4 , CHARS 
CHARACTER* 12 NAME 

COMMON/ BLKl /MESH 
COMMON/BLK9/CHAR2, CHAR3 , CHAR4 

C = 2. 997925E+08 
PI = 4. 0*ATAN(1. 0) 

READ(3,1070) NAME 
READ(3,1000) CHAR 
READ(3,1000) CHAR2 
READ(3,1000) CHARS 
READ(3,1000) CHAR4 
READ(3,1000) CHARS 
READ(3,*) MRESW 
READ(3,*) DISTW 
READ(3,*) DPER 
READ(3,*) LBIAS 
READ(3,*) GBIAS 
READ(3,*) MAXD 
READ(3,*) RES 
READ(3,*) NRES 
READ(3,*) METHOD 
READ(3,*) STARTND 
READ(3,*) STOPND 
READ(3,*) XORIGIN 
READ(3,*) YORIGIN 
READ(3,*) DD 
READ(3,*) A 
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c 



c 



c 



c 



c 



c 



READ(3,*) B 
ALPHA = CMPLX(A,B) 
READ(3,*) A 
READ(3,*) B 
BETA = CMPLX(A,B) 
READ(3,1000) CHARI 



IF( METHOD. EQ. 1) THEN 

WRITE(6,*) 'METHOD 1 SELECTED 
ELSEIF(METHOD. EQ. 2) THEN 

WRITE(6,*) 'METHOD 2 SELECTED 
ELSEIF(METHOD. EQ. 3) THEN 

WRITE(6,*) 'METHOD 3 SELECTED 
ELSE IF( METHOD. EQ. 4) THEN 

WRITE (6,*) 'METHOD 4 SELECTED 
ELSE IF( METHOD. EQ. 5) THEN 

WRITE(6,*) 'METHOD 5 SELECTED 

ELSE 

WRITE(6,*) 'METHOD 6 SELECTED 

ENDIF 



OBJECT BISECTION ' 

CONNECTED NODES ' 

SELECTED STOP NODE ' 
CONNECTED/SELECTED STOP NODE' 
EQUALLY SPACED CONNECTD NODE' 
SELECTED START AND STOP NODE' 



IF( ( CHAR2. EQ. ' I ' ) . OR. ( CHAR2. EQ. ' i ' ) ) THEN 

WRITE(6,*) 'INTERMEDIATE MATRIX FILE GENERATION - ENABLED' 

ELSE 

WRITE(6,*) 'INTERMEDIATE MATRIX FILE GENERATION - DISABLED' 

ENDIF 



IF((CHAR3. EQ. 'D').0R. (CHAR3.EQ. 'd')) THEN 

WRITE(6,*) 'DISSPLA FORTRAN PROGRAM GENERATION - ENABLED' 

ELSE 

WRITE(6,*) 'DISSPLA FORTRAN PROGRAM GENERATION - DISABLED' 

ENDIF 



IF((CHAR4.EQ. 'U').OR. (CHAR4. EQ. 'u')) THEN 

WRITE(6,*) 'UNIFORM MATERIAL SPECIFIED (NO INTERFACE) ' 

ELSE 

WRITE(6,*) 'MATERIAL SPECIFIED WITH A VACUUM AROUND OBJECT' 

ENDIF 



IF( (CHARS. EQ. 'M').0R. (CHARS. EQ. 'm')) THEN 

WRITE(6,*) 'MESH GENERATION «< ONLY »> - ENABLED' 

ELSE 

WRITE(6,*) 'mesh generation AND FE PROGRAM - ENABLED' 

ENDIF 

IF( (CHARI. EQ. 'P'). OR. (CHARI. EQ. 'p')) THEN 
READ(3,*) EO 
READ(3,*) FO 
LAMBDA = C/FO 
KO = 2*PI/ LAMBDA 

WRITE(6,*) 'PLANEWAVE BOUNDARY VALUE GENERATION - ENABLED' 
WRITE(6,*) 'AMPLITUDE(EO) = ', EO,', WAVENUMBER(KO) = ' ,K0 
WRITE( 6,*) 'WAVELENGTH = ', LAMBDA,', FREQUENCY(FO) = ' ,FO 
ELSEIF( (CHARI. EQ. 'C').OR. (CHARI. EQ. 'c')) THEN 
READ(3,*) EO 
RE AD (3,*) MODE 
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READ(3,*) FO 
LAMBDA = C/FO 
KO = 2*PI*LAMBDA 

WRITE(6,*) 'CLYINRICAL BOUNDARY VALUE GENERATION - ENABLED’ 
WRITE(6,*) ’AMPLITUDE(EO) = EO,', MODE NUMBER = MODE 

WRITE( 6,*) 'WAVELENGTH = ’.LAMBDA,', FREQUENCY(FO) = ’ ,FO 

ELSE 

WRITE(6,*) 'ALL BOUNDARY VALUE GENERATION METHODS - DISABLED’ 
DO 5, K = 1, RES 

READ(3,*) BCOND(K) 

5 CONTINUE 

END IF 



POLAR COORDINATE INPUT ROUTINE 

IF((CHAR. EQ. ’P’).OR. (CHAR. EQ. 'p')) THEN 
READ(3,1020) THETA, RADIUS 
SFAC = 2*PI*DD/RADIUS 

MESH(0,4) = (RADIUS*SIN(THETA*PI/180. 0))*SFAC + XORIGIN 
MESH(0,5) = (RADIUS*C0S(THETA*PI/180. 0))*SFAC + YORIGIN 
DO 10, J = 1, RES-1 

READ(3,1020) THETA, RADIUS 

MESH(J,4) = (RADIUS*SIN(THETA*PI/180. 0))*SFAC + XORIGIN 
MESH(J,5) = (RADIUS*COS(THETA*PI/180. 0))*SFAC + YORIGIN 
CONTINUE 

WRITE(6,*) 'POLAR COORDINATE INPUT SELECTED ’ 

RECTANGULAR COORDINATE INPUT ROUTINE 



20 



ELSEIF((CHAR. EQ. ’R’).0R. (CHAR. EQ. 'r')) THEN 
READ(3,1010) MESH(0,4), MESH(0,5) 

SFAC = 2*PI*DD/((MESH(0,4)**2+MESH(0,5)**2)**0. 5) 
MESH(0,4) = MESH(0,4)*SFAC + XORIGIN 

MESH(0,5) = MESH(0,5)*SFAC + YORIGIN 

DO 20, J = 1, RES-1 

READ(3,1010) MESH(J,4), MESH(J,5) 

MESH(J,4) = MESH(J,4)*SFAC + XORIGIN 

MESH(J,5) = MESH(J,5)*SFAC + YORIGIN 

CONTINUE 

WRITE(6,*) 'RECTANGULAR COORDINATE INPUT SELECTED ' 

ELSE 



WRITE(6,*) 'INPUT DATA FILE COORDINATE SPECIFICATION ERROR' 

END IF 

WRITE(*,1100) DD, SFAC 
MESH(RES,4) = MESH(0,4) 

MESH(RES,5) = MESH(0,5) 

NPNTS = J 

WRITE(6,*) 'NODE AT 0 DEGREES (X/Y COORDINATES) = ', MESH(0,4), 
CMESH(0,5) 

WRITE(6,*) 'X, Y OFFSETS = ', XORIGIN, YORIGIN 
WRITE(6,1090) ALPHA, BETA 
MRES = MRESW*2*PI 

WRITE(6,*) 'MESH RESOLUTION = ',MRESW,' WAVELENGTHS' 

DIST = DISTW*2*PI 

WRITE(6,*) 'CONTOUR DISTANCE = ' ,DISTW, ' WAVELENGTHS' 

WRITE(6,*) 'NUMBER OF DATA POINTS = ' ,RES 
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IF( (METHOD. EQ. 3). 
WRITE(6,*) ' 


OR. (METHOD. 


START NODE 


WRITE(6,*) ' 

ENDIF 


STOP NODE 


WRITE(4,*) '17' 


WRITE(4,*) '10 


1 


C 55 


1' 


WRITE(4,*) '2 


1 


C 80 


1' 


WRITE(4,*) '2 


1 


C 90 


2' 


WRITE(4,*) '2 


1 


C 100 


1' 


WRITE(4,*) '2 


1 


C 110 


2' 


WRITE(4,*) '2 


1 


C 120 


1' 


WRITE(4,*) '2 


1 


C 130 


2' 


WRITE(4,*) '2 


1 


C 140 


1' 


WRITE(4,*) '2 


1 


C 150 


2' 


WRITE(4,*) '2 


1 


C 160 


1' 


WRITE(4,*) '2 


1 


C 170 


2' 


WRITE(4,*) '2 


1 


C 180 


1' 


WRITE(4,*) '2 


1 


C 190 


2' 


WRITE(4,*) '2 


1 


C 200 


1' 


WRITE(4,*) '2 


1 


C 210 


2' 


WRITE(4,*) '2 


1 


C 220 


1' 


WRITE(4,*) '2 


1 


C 230 


2' 


WRITE(4,1070) NAME 
WRITE(4,*) 'L' 


WRITE (4,*) 'Mesh 
WRITE(4,*) 'L' 


Resolution 



WRITE(4,1030) MRESW 
WRITE(4,*) 'L' 

WRITE(4,*) 'Contour Distance' 
WRITE(4,*) 'L' 

WRITE(4,1030) DISTW 
WRITE(4,*) 'L' 

WRITE(4,*) 'Number of Points' 
WRITE(4,*) 'L' 

WRITE(4,1040) RES 
WRITE(4,*) 'L' 

WRITE(4,*) 'Method' 

WRITE(4,*) 'L' 



= ' , STARTND 
= ' , STOPND 



0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 



130 

365 

365 

365 

365 

365 

365 

365 

365 

365 

365 

365 

365 

365 

365 

365 

365 
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WRITE(4,1040) Method 
WRITE(4,*) 'L' 

WRITE(4,*) 'X Origin' 

WRITE(4,*) 'L' 

WRITE(4,1050) XORIGIN 
WRITE(4,*) 'L' 

WRITE(4,*) 'Y Origin' 

WRITE(4,*) 'L' 

WRITE(4,1050) YORIGIN 
WRITE(4,*) 'L' 

WRITE(4,*) 'Alpha' 

WRITE(4,*) 'L' 

WRITE(4,1060) ALPHA 
WRITE(4,*) 'L' 

WRITE(4,*) 'Beta' 

WRITE(4,*) 'L' 

WRITE(4,1060) BETA 
WRITE(4,*) 'L' 

WRITE(4,*) '9' 

C 

DO 30, II = 0, NPNTS-1 

WRITE(2,*) MESH(II,4), MESH(II,5) 

30 CONTINUE 

WRITE(2,*) MESH(0,4), MESH(0,5) 

WRITE(2,*) '999992, 999991 ' 

WRITE(2,*) '999990, 999990 ' 

C 

1000 FORMAT(Al) 

1010 F0RMAT(10X,F12. 8,4X,F12. 8) 

1020 FORMAT(3X,I3,5X,F12. 8) 

1030 FORMAT(F8. 6) 

1040 F0RMAT(I4) 

1050 F0RMAT(1X,F8. 6) 

1060 F0RMAT(1X,F8. 6, ' -J',F8.6) 

1070 FORMAT(A12) 

1080 FORMATC/) 

1090 FORMATC IX, 'ALPHA = ',F8.4,2X,'+ J ',F6. 4,', BETA = 'F8.4,2X,'+ J 
C ' ,F6. 4) 

1100 FORMATC IX, 'DESIRED DISTANCE = ',F8.5,', SCALE FACTOR = ',F8.5) 
RETURN 
END 
C 
C 

SUBROUTINE ROTATE C NPNTS , METHOD, STARTND,STOPND) 

C 

C THIS SUBROUTINE ROTATES THE SURFACE POINTS TO ALLOW FOR 
C A REARRANGING OF THE START NODE CFIRST DATA POINT). 

C 

REAL MESHCO: 200,5) 

INTEGER I, NPNTS, METHOD, STARTND, STOPND 
COMMON/BLKl/MESH 
C 
C 

IFCMETHOD. EQ. 6) THEN 

DO 10, I = 0, NPNTS-1 

MESHCI,!) = MESHCI,4) 



83 



o o o o o no 



MESH(I,2) = MESH(I,5) 

10 CONTINUE 

C 

DO 20, I = STARTND-1, NPNTS-1 

MESH(I-(STARTND-1),4) = MESH(I,1) 
HESH( I -(STARTND-1), 5) = MESH(I,2) 

20 CONTINUE 

C 

DO 30, 1=1, STARTND 

MESH(I+NPNTS-STARTND,4) = MESH( 1-1,1) 
MESHd+NPNTS- STARTND, 5) = MESH(I-1,2) 
30 CONTINUE 

END IF 

STOPND = STOPND - STARTND + 1 
IF(STOPND. GT. NPNTS) THEN 

STOPND = STOPND - NPNTS 
ELSEIFC STOPND. LT. 0) THEN 

STOPND = NPNTS + STOPND +1 
ELSEIFC STOPND. EQ. 0) THEN 
STOPND = 1 

ELSE 

STOPND = STOPND 

END IF 

RETURN 

END 



SUBROUTINE B OUND ( MRE S , NPNTS , METHOD , STOPND ) 

THIS SUBROUTINE REORDERS THE THE SURFACE POINTS 
BASED ON THE DESIRED INPUT MESH RESOLU INN. THE 
OBJECT IS ALSO BISECTED. 

REAL PERIM, DIST, MESH(0: 200,5) , TEMP, TEMPI, MRES, MRESN, MRESLW 
REAL MRESNL, MRESNR, TEMPR, TEMPL, DISTB, DZ, PI, MRESW, MRESRW 
INTEGER I, PERND, BIND, NPNTS , STARTND , STOPND, METHOD, NEWND, J 
INTEGER M, L 
COMMON/BLKl/MESH 
COMMON/BLK2/PERND,BIND 
C 

PI = 4. 0*ATAN(1. 0) 

PERIM = 0. 0 
C 

DO 10, I =0, NPNTS-1 

DIST = SQRT((MESH(I+1,4)-MESH(I,4))**2+(MESH(I+1,5) - 
C MESH(I,5))**2) 

PERIM = PERIM + DIST 
MESH( 1+1,3) = PERIM 
10 CONTINUE 

VRITE(6,*) 'PERIMETER LENGTH = ', PERIM 

PERND = NINT(PERIM/MRES - AMOD( PERIM/MRES , 2. 0) ) 

BIND = (PERND - 2)/2 

WRITE(6,*) 'PERIMETER NODE # = ', PERND 

MRESW = MRES/(2*PI) 

WRITE(6,*) 'YOUR REQUESTED MESH RESOLUTION OF ', MRESW 
IF((METHOD. EQ. 3). OR. (METHOD. EQ. 4).0R. (METHOD. EQ. 5). OR. 
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C (METHOD. EQ. 6)) THEN 

MRESNL = (PERIM -MESH( STOPND-1 ,3) )/FLOAT( PERND/2) 

MRESLW = MRESNL/ (2*PI) 

MRESNR = MESH( STOPND-1, 3)/FLOAT(PERND/2) 

MRESRW = MRESNR/ (2*PI) 

WRITE(6,*) . HAS BEEN MODIFIED TO . . . ' 

WRITE(6,*) 'LEFT SIDE OF SEGMENT . . . MRESLW 

WRITE(6,*) 'RIGHT SIDE OF SEGMENT . . . MRESRW 

ELSE 

MRESN = PERIM/FLOAT(PERND) 

MRESW = MRESN/(2*PI) 

WRITE(6,*) . HAS BEEN MODIFIED TO . . . MRESW 

ENDIF 

PERIMETER NODE INITIALIZATION (X,Y COORD) 

IF( (METHOD. EQ. 1). OR. (METHOD. EQ. 2)) THEN 
DO 30, 1=0, PERND-1 
TEMP = MRESN*I 
J = 1 

20 IF(TEMP.GT. MESH(J,3)) THEN 

J = J + 1 
GOTO 20 

ENDIF 

TEMPI = (TEMP - MESH(J-1,3))/(SQRT((MESH(J,4)-MESH(J-1,4)) 

C **2 + (MESH(J,5) - MESH(J-1,5))**2)) 

MESH(I+1,1) = MESH(J-1,4) + TEMP1*(MESH( J,4)-MESH( J-1,4)) 
MESH(I+1,2) = MESH(J-1,5) + TEMP1*(MESH( J,5)-MESH( J-1,5)) 

30 CONTINUE 

C 

ELSE 

DO 60, I = 0, PERND-1 
TEMPR = MRESNR*! 

TEMPL = PERIM - MRESNL*(PERND - I) 

IF(TEMPR. LE.MESH( STOPND-1, 3)) THEN 
J = 1 

40 IF(TEMPR. GT. MESH(J,3)) THEN 

J = J + 1 
GOTO 40 

ENDIF 

TEMPI = (TEMPR - MESH( J-1 , 3) )/( SQRT( (MESH( J,4) -MESH( J-1 ,4) ) 
C **2 + (MESH(J,5) - MESH(J-1,5))**2)) 

MESH(I+1,1) = MESH(J-1,4) + TEMP1*(MESH( J,4) -MESH( J-1 ,4) ) 
MESH(I+1,2) = MESH(J-1,5) + TEMP1*(MESH( J,5) -MESH( J-1 ,5) ) 
ELSE 

J = 1 

50 IF(TEMPL. GT. MESH(J,3)) THEN 

J = J + 1 
GOTO 50 

ENDIF 

TEMPI = (TEMPL - MESH( J-1, 3))/(SQRT((MESH(J,4)-MESH( J-1,4)) 
C **2 + (MESH(J,5) - MESH(J-1,5))**2)) 

MESH(I+1,1) = MESH(J-1,4) + TEMP1*(MESH(J,4)-MESH( J-1,4)) 
MESH(I+1,2) = MESH(J-1,5) + TEMP1*(MESH(J,5)-MESH( J-1,5)) 
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ENDIF 

60 CONTINUE 

ENDIF 
C 
C 

C BISECTION NODE INITIALIZATION (X,Y COORD) 

C 

WRITE(6,*) 'BISECTION NODE # = '.BIND 

IF( (METHOD. EQ. 1).0R. (METHOD. EQ. 3)) THEN 
DO 70, 1=1, BIND 

TEMPI = I/FLOAT(BIND + 1) 

MESH(I+PERND,1) = MESH(1,1) + TEMPl*(MESH(BIND+2 , 1) - 

CMESH(1,1)) 

MESH(I+PERND,2) = MESH(1,2) + TEMPl*(MESH(BIND+2,2) - 

CMESH(1,2)) 

70 CONTINUE 

C 

ELSE 

DO 80, 1=2, PERND+1 

MESH(PERND+I-1,1) = MESH(PERND+2-I , 1) + 0. 5*(MESH( I , 1) 
C MESH(PERND+2-I,l)) 

MESH(PERND+I-1,2) = MESH(PERND+2-I ,2) + 0. 5*(MESH( I ,2) 
C MESH(PERND+2-I,2)) 

80 CONTINUE 

C 

C 

ENDIF 

C 

IF( ( METHOD. EQ. 5 ) . OR. ( METHOD. EQ. 6 ) ) THEN 
DO 90, M = 1, BIND 

MESH(M,4) = MESH(PERND+M,1) 

MESH(M,5) = MESH(PERND+M,2) 

90 CONTINUE 

DISTB = 0.0 
MESH(0,3) =0.0 
DO 100, L = 1, BIND+1 
IF(L. EQ. 1) THEN 

DIST = SQRT((MESH(1,1) - MESH( 1 ,4) )**2 + 
C(MESH(1,2) - MESH(1,5))**2) 

ELSEIF(L. EQ. BIND+1) THEN 

DIST = SQRT((MESH(BIND,4) - MESH(BIND+2 , 1) ) 

C**2 + (MESH(BIND,5) - MESH(BIND+2 ,2) )**2) 

ELSE 

DIST = SQRT((MESH(L-1,4) - MESH(L,4)) 

C**2 + (MESH(L-1,5) - MESH(L,5))**2) 

ENDIF 

DISTB = DISTB + DIST 
MESH(L,3) = DISTB 
100 CONTINUE 

DZ = DISTB/(BIND + 1. 0) 

WRITE(6,*) 'DZ SPACING = ’ , DZ 

C 

DO 120, I = 1, BIND 
TEMP = DZ*I 
J = 1 

no IF(TEMP. GT. MESH(J,3)) THEN 
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j = j + 1 
GOTO 110 

ENDIF 

IF(J. EQ. 1) THEN 

TEMPI = TEMP/ (SQRT((MESH( 1,1) - MESH( 1 , 4) )**2 + 
C (MESH(1,2) - MESH(1,5))**2)) 

MESH(PERND+I,1) = MESH(1,1) + TEMP1*(MESH( 1 ,4) 

C - MESH(1,1)) 

MESH(PERND+I,2) = MESH(1,2) + TEMP1*(MESH( 1 ,5) 

C - MESH(1,2)) 

ELSE 

TEMPI = (TEMP - MESH(J-1,3))/(SQRT((MESH(J,4) - 

C MESH(J-1,4))**2 + (MESH(J,5) - MESH( J-1 , 5) )**2) ) 

MESH(PERND+I,1) = MESH(J-1,4) + TEMP1*(MESH( J,4) 
C - MESH(J-1,4)) 

MESH(PERND+I,2) = MESH(J-1,5) + TEMP1*(MESH( J,5) 
C - MESH(J-1,5)) 

ENDIF 

120 CONTINUE 

ENDIF 
RETURN 
END 



SUBROUTINE NORMAL 

THIS ROUTINE COMPUTES THE X AND Y COMPONENTS 
OF THE OUTWARD UNIT NORMAL AT EACH SURFACE POINT. 

REAL DR, DZ, DL, MESH( 0: 200 , 5) 

INTEGER I, PERND, BIND 
COMMON/BLKl/MESH 
C OMMON/ B LK2 / PERND , B I ND 
DO 10, 1=1, PERND 
IF(I. EQ. 1) THEN 

DR = MESH(2,1) - MESH( PERND, 1) 

DZ = MESH(2,2) - MESH(PERND,2) 

DL = SQRT(DR*DR+DZ*DZ) 

MESH(1,3)=-DZ/DL 
MESH(1,4)=DR/DL 
ELSEIFC I. EQ. PERND) THEN 

DR = MESH(1,1) - MESH( PERND- 1,1) 

DZ = MESH(1,2) - MESH( PERND- 1,2) 

DL = SQRT(DR*DR+DZ*DZ) 

MESH( PERND , 3)=-DZ/DL 
MESH(PERND,4)=DR/DL 

ELSE 

DR = MESH( 1+1,1) - MESH( 1-1,1) 

DZ = MESH(I+1,2)- MESH(I-1,2) 

DL = SQRT(DR*DR+DZ*DZ) 

MESH(I,3)=-DZ/DL 
MESH(I,4)=DR/DL 

ENDIF 
10 CONTINUE 
RETURN 
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END 



SUBROUTINE NODSET( METHOD , NABC , NBMX , UNK) 

THIS ROUTINE USES A TWO-SWEEP TECHNIQUE TO COMPUTE THE 
NUMBER OF NODES ALONG EACH NODE ROW, I. A MESH 
ORIENTATION ATTRIBUTE IS ALSO SET. 

SET Z-AXIS SPACING AND ENDPOINT NODES 

REAL DZ, ZZ, MESH(0; 200,5) , D, DIST, DISTB 
INTEGER I, J, NODES(200), MOR(200), PERND, NOLD, NNEW, BIND 
INTEGER L, METHOD, NABC( 100,3), NBMX, UNK 
CHARACTER*! CHAR2, CHAR3, CHAR4 
COMMON/BLKl/MESH 
COMMON/BLK2/PERND,BIND 
COMMON/BLK3/NODES 
COM.MON/BLK6/MOR 
COMMON/ BLK9/CHAR2 , CHAR3 , CHAR4 

UNK = 0 
NBMX = -99 

IF( (METHOD. EQ. l).OR. (METHOD. EQ. 3)) THEN 

DZ = (SQRT((MESH(1,2) - MESH(BIND+2,2))**2 + 

C (MESH(1,1) - MESH(BIND+2,1))**2))/(BIND + 1.0) 

ELSE 

DISTB =0.0 
DO 10, L = 1, BIND+1 
IF(L. EQ. 1) THEN 

DIST = SQRT((MESH(1,1) - MESH( PERND+1 , 1) )**2 + 
C(MESH(1,2) - MESH( PERND+1, 2 ))**2) 

ELSEIF(L. EQ. BIND+1) THEN 

DIST = SQRT((MESH(PERND+BIND,1) - MESH( BIND+2 , 1) ) 
C**2 + (MESH(PERND+BIND,2) - MESH( BIND+2 , 2) )**2) 

ELSE 

DIST = SQRT((MESH(PERND+I-1,1) - MESH( PERND+I , 1) ) 
C**2 + (MESH( PERND+I -1,2) - MESH( PERND+I ,2) )**2) 

ENDIF 

DISTB = DISTB + DIST 
10 CONTINUE 

DZ = DISTB/(BIND + 1. 0) 

ENDIF 

WRITE(6,*) 'BISECTION SPACING = ' , DZ 
NODES(l) = 2 
NODES(2) = 3 
NODES (BIND+2) = 2 
NODES (PERND+1) = 2 
NODES (PERND) = 3 
C PERFORMING FORWARD -SWEEP 
DO 20 1=3, BIND+1 

NOLD = NODES(I-l) 

D = SQRT((MESH(I,1) - MESH( I+PERND-1 , 1) )**2 + 

C ( MESH( 1,2) -MESH( I+PERND - 1 , 2 ) )**2 ) 

N.NEW = INT(0. 1 + D/DZ) + 2 
NODES(I) = NNEW 
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IF (NNEW. GT. NOLD+1) NODES(I) = NOLD + 1 
IF (NNEW. LT. NOLD- 1) NODES(I) = NOLD - 1 
IF (NODES(I). LT. 3) NODES(I) = 3 
20 CONTINUE 
C 

J = PER.ND + 2 

DO 30, I = PERND-1, BIND+3, -1 
NOLD = N0DES(I+1) 

D = SQRT((MESH(I,1) - MESH(J,1))**2 + 

C (MESH(I,2)-MESH(J,2))**2) 

NNEW = INT(0. 1 + D/DZ) + 2 
NODES(I) = NNEW 

IF (NNEW. GT. NOLD+1) NODES(I) = NOLD + 1 
IF (NNEW. LT. NOLD-1) NODES(I) = NOLD - 1 
IF (NODES(I). LT. 3) NODES(I) = 3 
J = J + 1 

30 CONTINUE 

C 

C BACK-SWEEP TO RESET LAST NODES IF NEEDED 

C 

I = BIND + 2 

40 1=1-1 

IF (I.EQ. 2) GO TO 50 

IF (NODES(I). LE. N0DES(I+1)+1) GO TO 60 
NODES(I) = NODES(I+l) + 1 
GO TO 40 
50 CONTINUE 

WRITE(6,*) ' PROGRAM ABORTED IN NODSET RIGHT SIDE BACKSWEEP ' 
STOP 

60 CONTINUE 

C 

I = BIND + 2 
70 1=1+1 

IF (I.EQ. PERND) GO TO 80 
IF (NODES(I). LE. N0DES(I-1)+1) GO TO 90 
NODES(I) = NODES(I-l) + 1 
GO TO 70 
80 CONTINUE 

WRITE(6,*) ' PROGRAM ABORTED IN NODSET LEFT HALF BACKSWEEP ' 
STOP 

90 CONTINUE 

C 

C FORWARD SWEEP TO LOAD MESH ORIENTATION ARRAY, MOR 

C 

MOR(l) = 0 
M0R(BIND+2) = 0 
M0R(PERND+1) = 0 
C 

DO 100, 1=2, BIND+1 

IF(N0DES(I+1).GT. NODES(D) THEN 
MOR(I) = 0 

ELSEIF(N0DES(I+1). LT.NODES(I)) THEN 
MOR(I) = 1 

ELSE 

IF(N0DES(I+2).GT. NODES(I)) THEN 
MOR(I) = 0 
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ELSEIF(N0DES(I+2). LT. NODES(I)) THEN 
MOR(I) = 1 

ELSE 

MOR(I) = MOR(I-l) 

ENDIF 

END IF 

100 CONTINUE 
C 

DO 110, I = PERND, BIND+3, -1 

IF(NODES(I-l).GT. NODES(I)) THEN 
MOR(I) = 0 

ELSEIF(NODES(I-l). LT. NODES(I)) THEN 
MOR(I) = 1 

ELSE 

IF(NODES(I-2).GT. NODES(D) THEN 
MOR(I) = 0 

ELSEIF(NODES(I-2). LT.NODES(I)) THEN 
MOR(I) = 1 

ELSE 

MOR(I) = MOR(I+l) 

ENDIF 

ENDIF 

110 CONTINUE 
C LOAD NABC ARRAY 

DO 120, I = 1, BIND+2 
IF(I.EQ. 1) THEN 

NABC(1,1) = 0 
NABC(1,2) = 1 
NABC(1,3) = 3 
ELSEIFCI. LE. BIND) THEN 

NABC(I,1) = NABC(I-1,2) 

NABC(I,2) = NABC(I-1,3) 

NABC(I,3) = NODES ( PERND -I+l) + N0DES(I+1) - 3 
ELSEIFCI.EQ. BIND+1) THEN 

NABC(I,1) = NABC(I-1,2) 

NABC(I,2) = NABC(I-1,3) 

NABC(I,3) = 1 

ELSE 

NABC(I,1) = NABC(I-1,2) 

NABC(I,2) = NABC(I-1,3) 

NABC(I,3) = 0 

ENDIF 

120 CONTINUE 
C 

DO 130, 1=1, BIND+2 

IF((CHAR2. EQ. ’ l' ). OR. (CHAR2. EQ. 'i')) THEN 

WRITE( 12,*) I ,NABC( I, 1) ,NABC( 1,2) ,NABC( 1,3) 

ENDIF 

UNK = UNK + NABC(I,2) 

IF(NABC(I,2). GE.NBMX) THEN 
NBMX = NABC(I,2) 

ENDIF 

130 CONTINUE 
C 

WRITE(*,*) 'MAXIMUM UNKNOWN WIDTH = ’ ,NBMX 
WRITEC*,*) 'TOTAL # OF UNKNOWNS = ' ,UNK 
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WRITE(*,*) 

RETURN 

END 



f I 



SUBROUTINE SORTER( I ,LEL,LAND) 

THIS SUBROUTINE GENERATES A MESH ROW FOR THE INPUT ROW I. 
LOADING OF THE LOCAL NODE-ELEMENT CONNECTION MATRICES LND AND 
NDL FOR ELEMENTS BETWEEN I AND I+l NODE ROWS. REFERENCE IS TO 
THE LEFT SIDE OF THE Ith ROW OR VECTOR. 



INTEGER N0DES(200), MOR(200), LND(0: 200,3) , NDL(200,4), NCT(200) 

INTEGER I, PERND, BIND, LEL, LAND, J, K, LL, JJ, NN, KK, N, Nl, N2 

INTEGER NDMX, NDSIL, NDS2L, NDSIR, NDS2R, LMX 

COMMON/BLK2/PERND,BIND 

COMMON/BLK3/NODES 

COMMON/ BLK4/ LND, NDL, NCT 

COMMON/BLK6/MOR 

IF(I.EQ. BIND+2) THEN 

WRITE (6,*) 'ERRORED OUT IN SUBROUTINE SORTER, YOU ATTEMPTED’ 
WRITE(6,*) ' TO CALL SORTER WITH I = BIND + 2! ' 

WRITE(6,*) ' THIS ROW HAS NO ELEMENTS' 

RETURN 



END IF 

DO 20, J = 0, 200 

DO 10, K = 1, 3 

LND(J,K) = 0 
CO.NTINUE 
CONTINUE 



NDMX = 200 

NDSIL = N0DES(PERND+2-I) 

NDS2L = NODES(PERND+l-I) 

NDSIR = NODES(I) 

NDS2R = NODESCI+1) 

LMX = NDSIL + NDSIR + NDS2L + NDS2R - 2 



LEFTSIDE OF THE BISECTION SEGMENT 
TOP ROW 

IFCI.EQ. 1) THEN 
LND(1,1) = 4 
LND(1,2) = 1 
LND(1,3) = 3 
LND(2,1) = 1 
LND(2,2) = 4 
LND(2,3) = 2 
LND(3,1) = 5 
LND(3,2) = 2 
LND(3,3) = 4 
LND(4,1) = 5 
LND(4,2) = 2 
LND(4,3) = 6 
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LND(5,1) 
LND(5,2) 
LND(5,3) 
LND(6,1) 
LND(6,2) 
LND(6,3) 
LEL = 6 
LAND = 7 



1 

6 

2 

6 

1 

7 



BOTTOM ROW 

ELSEIFd. EQ. BIND+1) THEN 
LND(1,1) = 2 
LND(1,2) = 6 
LND(1,3) = 1 
LND(2,1) = 6 
LND(2,2) = 2 
LND(2,3) = 7 
LND(3,1) = 3 
LND(3,2) = 7 
LND(3,3) = 2 
LND(4,1) = 3 
LND(4,2) = 7 
LND(4,3) = ^ 
LND(5,1) = 6 
LND(5,2) = 

L.ND(5,3) = 7 
LND(6,1) = ^ 
LND(6,2) = 6 
LND(6,3) = 5 
LEL = 6 
LAND = 7 

EQUAL NODE NUMBERS 



ELSEIFCNDSIL. EQ. NDS2L) THEN 
FOR MOR = 0 (LH ORIENTATION) 

IF(MOR(PERND+2-I).EQ. 0) THEN 



30 

40 

C 



FOR MOR = 
ELSE 



LND( 1,1) 
LND(1,2) 
LND(1,3) 
LND(2,1) 
LND(2,2) 
LND(2,3) 
DO 40, N 



= 1 



NDSIL + NDSIR 
2 

NDSIL + NDSIR 



+ 1 



2 

NDSIL + NDSIR 
1, NDSIL- 2 
N1 = 2*N + 1 
N2 = N1 + 1 
DO 30, K = 1, 3 

LND(N1,K) = LND(1,K) + N 
LND(N2,K) = LND(2,K) + N 
CONTINUE 
CONTINUE 

1 (RH ORIENTATION) 

LND(1,1) = NDSIL + NDSIR 
LND(1,2) = 1 
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LND(1,3) = NDSIL + NDSIR + 1 
LND(2,1) = 2 

LND(2,2) = NDSIL + NDSIR + 1 
LND(2,3) = 1 
DO 60, N = 1, NDSlL-2 
N1 = 2*N + 1 
N2 = N1 + 1 
DO 50, K = 1, 3 

LND(N1,K) = LND(1,K) + N 
LND(N2,K) = LND(2,K) + N 
50 CONTINUE 

60 CONTINUE 

ENDIF 
C 

C LEFTHAND MESH ORIENTATION 
C 

ELSEIF(MOR(PERND+2-I). EQ. 0) THEN 
LND(1,1) = NDSIL + NDSIR + 1 
LND(1,2) = 1 

LND(1,3) = NDSIL + NDSIR 
LND(0,1) = 0 

LND(0,2) = NDSIL + NDSIR 
LND(0,3) = 1 
DO 80, N = 1, NDSlL-1 
N1 = 2*N 
N2 = N1 + 1 
DO 70, K = 1, 3 

LND(N1,K) = LND(0,K) + N 
LND(N2,K) = LND(1,K) + N 
70 CONTINUE 

80 CONTINUE 

C 

C RIGHTHAND MESH ORIENTATION 

C 

ELSE 

LND(1,1) = 2 

LND(1,2) = NDSIL + NDSIR 
LND(1,3) = 1 

LND(0,1) = NDSIL + NDSIR - 1 
LND(0,2) = 1 

LND(0,3) = NDSIL + NDSIR 
DO 100, N = 1, NDSlL-2 
N1 = 2*N 
N2 = N1 + 1 
DO 90, K = 1, 3 

LND(N1,K) = LND(0,K) + N 
LND(N2,K) = LND(1,K) + N 
90 CONTINUE 

100 CONTINUE 
ENDIF 
C 

IF((I.EQ. l).OR. (I.EQ. BIND+1)) THEN 
GOTO 230 

ENDIF 

C 
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LEL = N2 



C 

C 

C 

C 

C 

C 

C 

C 



110 

120 



C 



130 

140 



C 

C 

C 



RIGHTSIDE OF THE BISECTION SEGMENT 



EQUAL NODE NUMBERS 



IF(NDS1R. EQ. NDS2R) THEN 
MOR = 0 (LH ORIENTATION) 

IF(MOR(I). EQ. 0) THEN 

LND(LEL+1,1) = NDSIL + NDSIR + NDS2L - 1 
LND(LEL+1,2) = NDSIL 
LND(LEL+1,3) = NDSIL + NDSIR + NDS2L 
LND(LEL+2,1) = NDSIL + 1 
LND(LEL+2,2) = NDSIL + NDSIR + NDS2L 
LND(LEL+2,3) = NDSIL 
DO 120, N = 1, NDSlR-2 
N1 = 2*N + 1 + LEL 
N2 = N1 + 1 

DO 110, K = 1, 3 

LND(N1,K) = LND(LEL+1,K) + N 
LND(N2,K) = LND(LEL+2,K) + N 
CONTINUE 
CONTINUE 

LEL = N2 
LAND = LMX 

MOR = 1 (RH ORIENTATION) 

ELSE 

LND(LEL+1,1) = NDSIL 

LND(LEL+1,2) = NDSIL + NDSIR + NDS2L - 1 
LND(LEL+1,3) = NDSIL + 1 
LND(LEL+2,1) = NDSIL + NDSIR + NDS2L 
LND(LEL+2,2) = NDSIL + 1 
LND(LEL+2,3) = NDSIL + NDSIR + NDS2L - 1 
DO 140, N = 1, NDSlR-2 
N1 = 2*N + 1 + LEL 
N2 = N1 + 1 

DO 130, K = 1, 3 

LND(N1,K) = LND(LEL+1,K) + N 
LND(N2,K) = LND(LEL+2,K) + N 
CONTINUE 
CONTINUE 

END IF \ 

LEL = N2 I 

LAND = LMX 



LEFT HAND MESH ORIENTATION 



ELSEIF(MORd). EQ. 0) THEN 
LND(LEL+1,1) = NDSIL 
LND(LEL+1,2) = NDSIL 
LND(LEL+1,3) = NDSIL 
LND(LEL+2,1) = NDSIL 
LND(LEL+2,2) = NDSIL 



+ 

+ 

+ 

+ 



NDSIR + NDS2L 

NDSIR + NDS2L 
1 

NDSIR + NDS2L 



1 



94 



c 



150 

160 

C 



170 

180 

C 



C 

C 

C 



C 



190 

200 

C 



210 

220 

C 



C 

C 

C 

C 

C 

C 

C 

230 



LND(LEL+2,3) = NDSIL 

DO 160, N = 1, NDSlR-1 
N1 = 2*N + LEL + 1 
DO 150, K = 1, 3 

LND(N1,K) = LND(LEL+1,K) + N 
CONTINUE 
CONTINUE 

DO 180, N = 1, NDSlR-2 
N2 = 2*N + LEL + 2 
DO 170, K = 1, 3 

LND(N2,K) = LND(LEL+2,K) + N 
CONTINUE 
CONTINUE 

LEL = N1 
LAND = LMX 



RIGHT HAND MESH ORIENTATION 



ELSE 

LND(LEL+1,1) = 
LND(LEL+1,2) = 
LND(LEL+1,3) = 
LND(LEL+2,1) = 
LND(LEL+2,2) = 
LND(LEL+2,3) = 



NDSIL 

NDSIL + NDSIR + NDS2L 
NDSIL + 1 

NDSIL + NDSIR + NDS2L 
NDSIL + 1 

NDSIL + NDSIR + NDS2L 



1 



1 



DO 200, N = 1, NDSlR-2 
N1 = 2*N + LEL + 1 
DO 190, K = 1, 3 

LND(N1,K) = LND(LEL+1,K) + N 
CONTINUE 
CONTINUE 



DO 220, N = 1, NDSlR-3 
N2 = 2*N + LEL + 2 
DO 210, K = 1, 3 

LND(N2,K) = LND(LEL+2,K) + N 
CONTINUE 
CONTINUE 

LEL = N1 
LAND = LMX 



ENDIF 



ZEROING NDL AND NCT PRIOR TO FILL 

DO 250, N = 1, NDMX 
NCT(N) = 0 
DO 240, K = 1, 4 
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240 

250 

C 

C 

C 



260 

270 

C 

C 

C 

C 

C 

C 

C 



NDL(N,K) = 0 
CONTINUE 
CONTINUE 

SCANNING END ARRAY TO LOAD NDL 

DO 270 LL = 1, LEL 

DO 260 JJ = 1, 3 

NN = LND(LL,JJ) 
NCT(NN) = NCT(NN) + 1 
KK = NCT(NN) 
NDL(NN,KK) = LL 
CONTINUE 
CONTINUE 

END 



SUBROUTINE FINDER( I , LAND ,DIST) 

THIS SUBROUTINE DETERMINES THE (X,Y) COORDINATES OF EACH 
NODE IN THE CALLING Ith ROW OR VECTOR MESH. 



INTEGER I, J, LAND, PERND, BIND, NODES(200) 

REAL MESH(0: 200,5) , NDP(200,2), DIST 

COMMON/BLKl/MESH 

COMMON/BLK2/PERND , BIND 

COMMON/BLK3/NODES 

COMMON/BLK5/NDP 

IF(I. EQ. 1) THEN 

NDP(1,1) = MESH(1,1)+MESH(1,3)*DIST 
NDP(1,2) = MESH(1,2)+MESH(1,4)*DIST 
NDP(2,1) = MESH(1,1) 

NDP(2,2) = MESH(1,2) 

NDP(3,1) = MESH( PERND, 1)+MESH( PERND, 3)*DIST 
NDP(3,2) = MESH(PERND,2)+MESH(PERND,4)*DIST 
NDP(4,1) = MESH( PERND, 1) 

NDP(4,2) = MESH(PERND,2) 

NDP(5,1) = MESH(PERND+1,1) 

NDP(5,2) = MESH(PERND+1,2) 

NDP(6,1) = MESH(2,1) 

NDP(6,2) = MESH(2,2) 

NDP(7,1) = MESH(2,10+MESH(2,3)*DIST 
NDP(7,2) = MESH(2,2)+MESH(2,4)*DIST 



ELSEIFd. EQ. BIND+1) THEN 

NDP(1,1) = MESH(BIND+3,1)+MESH(BIND+3,3)*DIST 
NDP(1,2) = MESH(BIND+3,2)+MESH(BIND+3,4)*DIST 
NDP(2,1) = MESH(BIND+3,1) 

NDP(2,2) = MESH(BIND+3,2) 

NDP(3,1) = MESH(PERND+BIND,1) 

NDP(3,2) = MESH(PERND+BIND,2) 

NDP(4,1) = MESH( BIND+1, 1) 

NDP(4,2) = MESH( BIND+1, 2) 

NDP(5,1) = MESH( BIND+1, 1)+MESH( BIND+1, 3)*DIST 
NDP(5,2) = MESH( BIND+1, 2)+MESH( BIND+1, 4)*DIST 
NDP(6,1) = MESH(BIND+2, 1)+MESH(BIND+2,3)*DIST 
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NDP(6,2) = MESH(BIND+2,2)+MESH(BIND+2,4)*DIST 
NDP(7,1) = MESH(BIND+2,1) 

NDP(7,2) = MESH(BIND+2,2) 

ELSEIFCI. EQ. BIND+2) THEN 

WRITE(6,*) ' ERRORED OUT IN SUBROUTINE FINDER, YOU ATTEMPTED' 
WRITE(6,*) ' TO CALL FINDER WITH I = BIND + 2! ' 

WRITE(6,*) ' THIS ROW HAS NO ELEMENTS AND NO COORDINATES' 

ELSE 

C NODE 1 

NDP(I,I) = MESH(PERND-I+2,I)+MESH(PERND-I+2,3)*DIST 
NDP(1,2) = MESH(PERND-I+2,2)+MESH(PERND-I+2,4)*DIST 
C NODE 2 TO THE BISECTION SEGMENT 

DO 10, J = 2, NODE SCPERND- 1+2) 

NDP(J,1)=MESH(PERND-I+2,1)-(J-2)*(MESH(PERND-I+2,1)- 
C ME SH ( PERND+ 1 - 1 , 1 ) ) / ( NODES ( PERND - 1 +2 ) - 2 ) 

NDP(J,2)=MESH(PERND-I+2,2)-(J-2)*(MESH(PERND-I+2,2)- 
C MESH(PERND+I-l,2))/(N0DES(PERND-I+2)-2) 

10 CONTINUE 

C BISECTION SEGMENT TO THE RIGHTSIDE SURFACE 

DO 20, J = 3, NODES(I) 

NDP(NODES(PERND-I+2)+J-2,l)=MESH(PERND+I-l,l)+(J-2)* 

C ( MESH( 1,1) -MESH( PERND+I -1,1) ) /( NODES( I ) -2) 

NDP ( NODE S ( PERND - 1 +2 ) + J - 2 , 2 ) =ME S H ( PE RND+ I-l,2)+(J-2)* 

C ( MESH( I , 2 ) -MESH( PERND+I - 1 , 2 ) ) / ( NODES( I ) - 2 ) 

20 CONTINUE 

C Ith ROWS LAST NODE 

NDP(NODES(PERND-I+2)+NODES( I) -1 , 1)=MESH( I , 1)+MESH( I ,3)*DIST 
NDP(NODES(PERND-I+2)+NODES( I) -1 ,2)=MESH( I ,2)+MESH( I ,4)*DIST 
C I+lth ROWS FIRST NODE 

NDP( NODES( PERND- 1+2 )+NODES( I ) , 1 )=MESH( PERND- I+l , 1 )+MESH 
C (PERND-I+1,3)*DIST 

NDP(NODES(PERND-I+2)+NODES(I),2)=MESH(PERND-I+l,2)+MESH 
C (PERND-I+1,4)*DIST 

C I+ITH ROW (NODE 2) TO THE BISECTION SEGMENT 

DO 30, J = 2, NODES(PERND-I+l) 

NDP( J+NODESC PERND-I+2)+NODES( I ) - 1 , 1)=MESH( PERND- I+l , 

C 1)-(J-2)*(MESH(PERND-I+1,1)-MESH( PERND+I, 1))/ 

C (NODES(PERND-I+l)-2) 

NDP ( J+NODE S ( PERND - 1 +2 ) +NODE S ( I ) - 1 , 2 ) =ME SH( PERND - 1 +1 , 

C 2)-(J-2)*(MESH(PERND-I+l,2)-MESH(PERND+I,2))/ 

C (NODES(PERND-I+l)-2) 

30 CONTINUE 

C I+lth ROW BISECTION SEGMENT TO THE RIGHTSIDE SURFACE 

DO 40, J = 3, NODES(I+l) 

NDP(LAND-NODES(I+l)+J-l,l)=MESH( PERND+I, l)+(J-2)* 

C (MESH( I+l , 1) -MESH( PERND+I , 1) )/(NODES( I+l) -2) 

NDP( LAND-NODESC I+l)+J- 1 , 2)=MESH( PERND+I , 2)+( J-2)* 

C ( MESH( 1+1,2) -MESH( PERND+I , 2 ) ) / ( NODES( 1+ 1 ) - 2 ) 

40 CONTINUE 

C LAST NODE 

NDP(LAND,1) = MESH(I+1,1)+MESH(I+1,3)*DIST 
NDP(LAND,2) = MESH( I+l ,2)+MESH( I+l ,4)*DIST 
C 

ENDIF 

RETURN 
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END 



SUBROUTINE VARINT( J , F , ALPHA , BETA , AREA , LND) 

GENERATING VARIATIONAL FINITE ELEMENT AREA INTEGRATIONS OF THE 
LINEAR BASIS FUNCTION LAGRANGIAN FOR THE HELMHOLTZ EQUATION. 
THESE ARE RETURNED IN F(3,3). X(3) AND Y(3) ARE THE WAVENUMBER 
NORMALIZED CARTESIAN COORDINATES OF THE TRIANGLE VERTICES. 

X = Ko*x, Y = Ko*y ALPHA AND BETA ARE COMPLEX 
MATERIAL PARAMETERS WITHIN THE ELEMENT. 

FOR TM INCIDENCE: ALPHA = 1/ur; BETA = er 
FOR TE INCIDENCE: ALPHA = 1/er; BETA = ur 

OUTPUTS : F(3,3) - FINITE ELEMENT AREA INTEGRATION 

AREA - AREA OFF A TRIANGLE 

COMPLEX ALPHA, BETA, B12, F(3,3) 

REAL NDP( 200,2) , X(3), Y(3), T(3,3), AREA, DET 
INTEGER L, K, J, LND(0: 200,3) 

COMMON/ BLK5/NDP 

X(l) = NDP(LND(J,1),1) 

X(2) = NDP(LND(J,2),1) 

X(3) = NDP(LND(J,3) , 1) 

Y(l) = NDP(LND(J,1) ,2) 

Y(2) = NDP(LND(J,2) ,2) 

Y(3) = NDP(LND(J,3) ,2) 

DET = X(2)*Y(3) + X(3)*Y(1) + X(1)*Y(2) - X(3)*Y(2) - 
CX(1)*Y(3) - X(2)*Y(1) 

AREA = ABS(0.5*DET) 

B12 = BETA/12. 

T(l,l) = (Y(2) - Y(3))/DET 

T(l,2) = (Y(3) - Y(1))/DET 

T(l,3) = (Y(l) - Y(2))/DET 

T(2,l) = (X(3) - X(2))/DET 

T(2,2) = (X(l) - X(3))/DET 

T(2,3) = (X(2) - X(1))/DET 

T(3,l) = (X(2)*Y(3) - X(3)*Y(2))/DET 

T(3,2) = (X(3)*Y(1) - X(1)*Y(3))/DET 

T(3,3) = (X(1)*Y(2) - X(2)*Y(1))/DET 

DO 10, K = 1, 3 

DO 10, L = 1, 3 

F(K,L) = ALPHA*(T(1,K)*T(1,L) + T( 2 , K)*T( 2 , L) ) - B12 
IF(K.EQ.L) F(K,L) = F(K,L) - B12 
F(K,L) = AREA*F(K,L) 

RETURN 

END 



SUBROUTINE L0ADER( BCOND , OFFSET , ALPHA , BETA , NABC , IMX , NBMX , EO , SURBC , 
CCHARl , LINE , MODE , XORIGIN , YORIGIN , SVEC , CHARS ) 

C 

COMPLEX A(50,50),B(50,50),C(50,50),P(50,100) 
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COMPLEX F(3,3),FR0W(100,3,3) ,BCOND(100) ,LINE(50) 

COMPLEX ALPHA, BETA, DETERM, SURBC(IOO), SVEC(50,100) 

REAL OFFSET , MINAREA , MAXAREA , AREA , RATIO , EO , KO , XORIGIN , YORIGIN 
REAL UBCOND 

INTEGER I , J , JD , K , L , NDTOP , NDBOT , NDTOT , NOD , KND , ND( 3 ) , LEL , LAND 
INTEGER LND(0: 200,3) , NDL(200,4), NCT(200), PERND, BIND, JJ 
INTEGER NODES(200), MINROW, MINEL, MAXROW, MAXEL, TCALL, NL 
INTEGER N,M,NBMX,NABC( 100,3) ,INORM,IMX, MODE 
CHARACTER*! CHARI, CHAR2, CHAR3, CHAR4, CHARS 
C 

C0MM0N/BLK2/PERND , BIND 

C0MM0N/BLK3/N0DES 

C0MM0N/BLK4/LND,NDL,NCT 

C0MM0N/BLK7 /MINAREA , MINROW , MINEL , MAXAREA , MAXROW , MAXEL , AREA 
COMMON/BLK8/A,B,C,P 
COMMON/BLK9/CHAR2, CHAR3, CHAR4 
C 

UBCOND =1.0 
I NORM = 0 
IMX = BIND + 2 
I = 1 

TCALL = BIND + 1 
C 

IF((CHAR3.EQ. 'D’).0R. (CHAR3.EQ. 'd')) THEN 
WRITE(20,1030) 

WRITE(20,1040) 

WRITE(20,1050) 

WRITE(20, 1060) 

WRITE(20,1070) 

ENDIF 

WRITE(*,1000) I, TCALL 
C 

CALL SORTER(I,LEL,LAND) 

CALL FINDERC I, LAND, OFFSET) 

IF(. NOT. ((CHARS. EQ. ’M’ ). OR. (CHARS. EQ. 'tn'))) THEN 
CALL BNDC( I , BCOND , EO , SURBC , CHARI , ALPHA , BETA , LINE , MODE , XORIGIN , 
CY0RIGIN,CHAR4) 

ENDIF 

CALL FILL( I , LEL , FROW , ALPHA , BETA ) 

IF(.NOT. ((CHARS. EQ. 'M').0R. (CHARS. EQ. 'm'))) THEN 

CALL ZERO 

ENDIF 

C ESTABLISH THE NUMBER OF NODES IN THE I = 1 AND I = 2 ROWS 
NDTOP = 2 
NDBOT = S 
NDTOT = 7 

B, C & P FOR I = 1 (TOP) 

J = 1 

DO 70, JD = 1, NCT(2) 

L = NDL(2,JD) 

DO SO, K = 1, 3 

ND(K) = LND(L,K) 

IF(ND(K).EQ. 2) KND = K 
SO CONTINUE 
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DO 60, K = 1, 3 
NL = ND(K) 

C BOUNDARY CONDITION 

IF(NL. EQ. 1) THEN 

P(J,1) = -UBCOND*FROW(L,KND,K) + P(J,1) 

C UPPER ROW 

ELSEIF(NL. EQ. 2) THEN 

B(J,1) = FROW(L,KND,K) + B(J,1) 

C LOWER ROW 

ELSEIF((NL. GT. 3). AND. (NL. LT. 7)) THEN 

C(J,NL-3) = FROW(L,KND,K) + C(J,NL-3) 

C ERROR 

ELSE 

WRITE(*,*) NL, 'ERROR - INDEX EQUALS 3 OR 7' 

ENDIF 

60 CONTINUE 

70 CONTINUE 

C 

IF((CHAR2. EQ. ' I ' ) . OR. ( CHAR2. EQ. 'i')) THEN 
WRITEd,*) I 
DO 72, M = 1, NABC(I,2) 

WRITE(1,1020) (REAL(B(M,N)) , N= 1, NABC(I,2)) 

72 CONTINUE 

WRITEd,*) ' ' 

DO 73, M = 1, NABC(I,2) 

WRITEd, 1020) (REAL(C(M,N)), N = 1, NABC(I,3)) 

73 CONTINUE 

WRITEd,*) ’ ' 

DO 74, M = 1, NABC(I,2) 

WRITE(1,1020) (REAL(P(M,N)) , N = 1, PERND) 

74 CONTINUE 
C 

WRITEd,*) ' ' 

WRITE(1,*) ' ' 

WRITEd,*) ’ ' 

ENDIF 

C 

IF(. NOT. ((CHARS. EQ. 'M' ). OR. (CHARS. EQ. 'm'))) THEN 

CALL MARCH( I , IMX , NBMX , NABC ( I , 1 ) , NABC ( I , 2 ) , NABC (1,3), INORM , S VEC ) 

CALL ZERO 

ENDIF 

C BEGIN LOOPING 

DO 900, 1=1, BIND 
C 

DO 400, J = 1, NDBOT-2 

NOD = J + NDTOP + 1 
DO 300, JD = 1, NCT(NOD) 

L = NDL(NOD,JD) 

DO 100, K = 1, 3 

ND(K) = LND(L,K) 

IF(ND(K).EQ. NOD) KND = K 
100 CONTINUE 

DO 200, K = 1, 3 
NL = ND(K) 

C BOUNDARY CONDITION 

IF((I.EQ. 1). AND. (NL. EQ. D) THEN 
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P(J,NL) = -UBCOND*FROW(L,KND,K) + P(J,NL) 

ELSEIFCNL. EQ. 1) THEN 

P ( J , PERND - 1+2 ) = -UBCOND*FROW ( L , KND , K ) +P ( J , PERND - I +2 ) 

C SPECIAL CASE FOR NODE j'/2 AND I = 1 
ELSEIF((I.EQ. D.AND. (NL. EQ. 2))THEN 
A(J,1) = FROW(L,KND,K) + A(J,1) 

ELSEIF(NL. EQ. NDTOP) THEN 

P(J,I) = -UBCOND*FROW(L,KND,K) + P(J,I) 

ELSEIFCNL. EQ. NDTOP+1) THEN 

P(J,PERND-I+l)=-UBCOND*FROW(L,KND,K)+P(J,PERND-I+l) 
ELSEIFCNL. EQ.NDTOT) THEN 

PCJ,I+1) = -UBCOND*FROWCL,KND,K) + PCJ,I+1) 

C UPPER ROW 

ELSEIFCNL. LT. NDTOP) THEN 

ACJ,NL-1) = FROWCL,KND,K) + ACJ,NL-1) 

C LOWER ROW 

ELSEIFCNL. LT. NDTOT) THEN 

BCJ,NL-NDTOP-l) = FROWCL,KND,K) + BC J,NL-NDTOP-l) 

C ERROR 

ELSE 

WRITEC*,*) NL, 'ERROR - DUE TO INDEX GREATER THAN NODE TOTAL’ 

END IF 
C 

200 CONTINUE 

300 CONTINUE 

400 CONTINUE 

INDEX FOR NEXT ROW AND COMPLETE B, C, P FILLS 

JJ = I + 1 

WRITEC6,1010) JJ,TCALL 
CALL SORTERCJJ,LEL,LAND) 

CALL FINDERCJJ,LAND,OFFSET) 

IFC. NOT. CC CHARS. EQ. 'M’).0R. C CHARS. EQ. 'm'))) THEN 
CALL BNDCC J J , BCOND , EO , SURBC , CHARI , ALPHA , BETA , LINE , MODE , 

C X0RIGIN,Y0RIGIN,CHAR4) 

END IF 

CALL F I LLC J J , LE L , FROW , ALPHA , BETA ) 

C ESTABLISH THE NUMBER OF NODES IN THE I+lth AND I+2th ROWS 
IFCJJ. NE. CBIND+1)) THEN 

NDTOP = N0DESCPERND+2-JJ) + NODESCJJ) - 1 
NDBOT = NODESCPERND+l-JJ) + NODESCJJ+1) - 1 
NDTOT = NDBOT + NDTOP 

ELSE 

NDTOP = S 
NDBOT = 2 
NDTOT = 7 

ENDIF 

C 

DO 800, J = 1, NDTOP-2 
NOD = J + 1 

DO 700, JD = 1, NCTCNOD) 

L = NDLCNOD,JD) 

DO SOO, K = 1, 3 

NDCK) = LNDCL,K) 
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IF(ND(K).EQ. NOD) KND = K 
500 CONTINUE 

DO 600, K = 1, 3 
NL = ND(K) 

C BOUNDARY CONDITION 
IF(NL. EO-1) THEN 

P ( J , PERND - J J+2 ) = - UBC0ND*FR0W( L , KND , K ) +P ( J , PERND - J J+2 ) 
ELSEIF(NL. EO. NDTOP) THEN 

P(J,JJ) = -UBCOND*FROW(L,KND,K) + P(J,JJ) 

ELSEIF(NL. EO- NDTOP+1) THEN 

P ( J , PERND - J J+ 1 ) = - UB COND*FROW( L , KND , K ) +P ( J , PERND - J J+ 1 ) 
ELSEIF(( JJ. EO. (BIND+1)). AND. (NL. EO. NDTOT)) THEN 
C(J,1) = FROW(L,KND,K) + C(J,1) 

ELSEIF(NL.EO. NDTOT) THEN 

P(J,JJ+1) = -UBCOND*FROW(L,KND,K) + P(J,JJ+1) 

C UPPER ROW 

ELSEIF(NL. LT. NDTOP) THEN 

B(J,NL-1) = FROW(L,KND,K) + B(J,NL-1) 

C LOWER ROW 

ELSEIF(NL. LT. NDTOT) THEN 

C(J,NL-NDT0P-1) = FROW(L,KND,K) + C( J,NL-NDT0P-1) 

C ERROR 

ELSE 

WRITE(*,*) NL, 'ERROR - DUE TO INDEX GREATER THAN NODE TOTAL 
ENDIF 
C 

600 CONTINUE 

700 CONTINUE 

800 CONTINUE 

C 

IF( ( CHAR2. EO. ' I ’ ) . OR. ( CHAR2. EO. ' i ' ) ) THEN 
WRITEd,*) JJ 
DO 91, M = 1, NABC(JJ,2) 

WRITEd, 1020) (REAL(A(M,N)), N = 1, NABC(JJ,1)) 

91 CONTINUE 

WRITEd,*) ' ’ 

DO 92, M = 1, NABC(JJ,2) 

WRITEd, 1020) (REAL(B(M,N)) , N = 1, NABC(JJ,2)) 

92 CONTINUE 

WRITEd,*) ' ’ 

DO 93, M = 1, NABC(JJ,2) 

WRITEd, 1020) (REAL(C(M,N)), N = 1, NABC(JJ,3)) 

93 CONTINUE 

WRITE(1,*) ' ' 

DO 94, M = 1, NABC(JJ,2) 

WRITEd, 1020) (REAL(P(M,N)), N = 1, PERND) 

94 CONTINUE 
C 

WRITEd,*) ' ' 

WRITEd,*) ’ ’ 

WRITEd,*) ' ' 

ENDIF 

C 

IF(. NOT. (( CHARS. EO. ' M' ). OR. (CHARS. EO- 'm'))) THEN 
CALL MARCH( JJ , IMX , NBMX , NABC ( JJ , 1 ) , NABC( JJ , 2 ) , NABC( JJ , 3 ) , INORM 
CSVEC) 
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CALL ZERO 
ENDIF 
C 

900 CONTINUE 

LOAD A, B & P FOR I = BIND + 2 (BOTTOM) 

J = 1 

DO 30, JD = 1, NCT(7) 

L = NT)L(7,JD) 

DO 10, K = 1, 3 

ND(K) = LND(L,K) 

IF(ND(K).EQ. 7) KND = K 
10 CONTINUE 

DO 20, K = 1, 3 
NL = ND(K) 

C BOUNDARY CONDITION 

IF(NL.EQ. 6) THEN 

P(J,JJ+1) = -UBCOND*FROW(L,KND,K) + P(J,JJ+1) 

C UPPER ROW 

ELSEIF((NL.GT. 1). AND. (NL. LT. 5)) THEN 

A(J,NL-1) = FROW(L,KND,K) + A(J,NL-1) 

C LOWER ROW 

ELSEIF(NL.EQ. 7) THEN 

B(J,1) = FROW(L,KND,K) + B(J,1) 

C ERROR 

ELSE 

WRITE(*,*) NL, 'ERROR - INDEX EQUAL TO 1 OR 5’ 

ENDIF 

C 

20 CONTINUE 

30 CONTINUE 

C 

IF((CHAR2. EQ. 'l').OR. (CHAR2.EQ. 'i')) THEN 
WRITE (1,*) TCALL+1 
DO 61, M = 1, NABC(TCALL+1,2) 

WRITE( 1,1020) (REAL(A(M,N)), N = 1, NABC(TCALL+1 , 1) ) 

61 CONTINUE 

WRITEd,*) ' ’ 

DO 62, M = 1, NABC(TCALL+1,2) 

WRITE(1,1020) (REAL(B(M,N)) , N = 1, NABC(TCALL+1 , 2) ) 

62 CONTINUE 

WRITEd,*) ’ ' 

DO 64, M = 1, NABC( TCALL+1, 2) 

WRITEd, 1020) (REAL(P(M,N)), N = 1, PERND) 

64 CONTINUE 

C 

WRITEd,*) ' ’ 

WRITEd,*) ' ' 

WRITEd,*) ' ’ 

ENDIF 

C 

WRITE(6,*) ' FINAL INVERSION' 

IF(. NOT. ((CHARS. EQ. ' M' ). OR. (CHARS. EQ. 'm'))) THEN 
CALL MARCH( J J+ 1 , IMX , NBMX , NABC( J J+ 1 , 1 ) , NABC ( J J+ 1 , 2 ) , NABC ( J J+ 1 , 3 ) , 
CINORM,SVEC) 
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c 

c 

1000 

1010 

1020 

1030 

1040 

1050 

1060 

1070 

1080 

1090 

1100 

C 

C 

C 

C 

C 



END IF 



,M INAREA 

AND ELEMENT NUMBER 



ELEMENT NUMBER 



WRITE(2,*) ’ END END 
WRITE(6,*) 

WRITE (6,*) 

WRITE(6,*) 

WRITE (6,*) 

WRITE(6,*) 

RATIO = MAXAREA/M INAREA 
WRITE(6,*) 'AREA RATIO = '.RATIO 
IFCRATIO. GT. 2. 5) THEN 

WRITE(6,*) 'YOU SHOULD CONSIDER ABORTING THIS RUN AND ' 
WRITE(6,*) 'LOOKING AT THE MESH IN CURVE DIGITIZER ' 
WRITE(6,*) 'A BETTER METHOD MAY BE AVAILABLE ' 



MINIMUM AREA = 

AT ROW ' .MINROW, ' 

MAXIMUM AREA = ‘ .MAXAREA 
AT ROW ' .MAXROW, ' AND 



' ,MINEL 
' .MAXEL 



ENDIF 

IF( ( CHARS. EQ. ' D ' ) . OR. ( CHARS. EQ. ' d ' ) ) THEN 
WRITE(20,1080) 

WRITE(20,1090) 

WRITE(20, 1100) 

ENDIF 



FORMATC 16,' OUT OF', 16,' CALLS ') 
FORMATC 16,' OUT OF', 16,' CALLS ') 
FORMATC 20CE14. 8,1X,E14. 8,1X)) 
FORMATC 6X, ' CALL COMPRS ' ) 

FORMATC 6X, 'CALL NOBRDR' ) 

FORMATC 6X,'CALL PAGEC 8. 0 , 10. 0) ' ) 
FORMATC 6X,'CALL AREA2DC5. 0 , 7. 0) ' ) 
FORMATC 6X,'CALL FRAME') 

FORMATC 6X,'CALL DONEPL' ) 

FORMATC 6X,'STOP') 

FORMATC 6X,'END') 



RETURN 

END 



SUBROUTINE FILLC I , LEL , FROW , ALPHA , BETA) 

COMPLEX FC3,S), FROWC 100 , 3 , 3) , ALPHA, BETA, A, B 
REAL AREA, MINAREA, MAXAREA 
REAL NDPC200,2) 

INTEGER I, J, K, L, LEL, LNDC 0: 200 , 3) , NDLC200,4), NCTC200) 

INTEGER MINROW, MINEL, MAXROW, MAXEL, M 

CHARACTER*! CHAR2, CHARS, CHAR4 

C0MM0N/BLK4/LND,NDL,NCT 

COMMON/BLK5/NDP 

C0MM0N/BLK7 /MINAREA , MINROW , MINEL , MAXAREA , MAXROW , MAXEL , AREA 
COMMON/ BLK9/CHAR2, CHARS, CHAR4 

IFC C CHAR2. EQ. ' I ' ) . OR. C CHAR2. EQ. ' i ' ) ) THEN 
WRITE C 11,*) 'ROW NUMBER = ',I 

ENDIF 

DO 30, J = 1, LEL 

IFCCCHAR4.EQ. 'U').OR. CCHAR4. EQ. 'u')) THEN 
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10 



20 



C 



C 

30 

C 

1000 

1010 

1020 

C 



A = ALPHA 
B = BETA 

ELSE 

IF(( J. EQ. 1). OR. ( J. EQ. 2). OR. (J. EQ. LEL-1). OR. ( J. EQ. LEL)) 
C THEN 

A = (1. 0,0. 0) 

B = ( 1. 0,0. 0) 

ELSE 

A = ALPHA 
B = BETA 

END IF 



END IF 

CALL VARINT(J,F,A,B,AREA,LND) 

IF( ( J. GT. 2 ) . AND. ( J. LT. LEL- 1 ) ) THEN 

IF(AREA. LT. MINAREA) THEN 
MINAREA = AREA 
MINROW = I 
MINEL = J 

ELSEIFCAREA. GT. MAXAREA) THEN 
MAXAREA = AREA 
MAXROW = I 
MAXEL = J 

END IF 



ENDIF 



DO 20, K = 1, 3 

WRITE(2,*) NDP(LND(J,K),1), NDP(LND( J,K) ,2) 

DO 10, L = 1, 3 

FROW(J,K,L) = F(K,L) 

CONTINUE 

IF((CHAR2. EQ. ' I ’ ). OR. (CHAR2. EQ. 'i')) THEN 

WRITE(11,1000) J,K,(REAL(F(K,L)), L = 1,3) 

ENDIF 

CONTINUE 

IF((CHAR2.EQ. 'l’).OR. (CHAR2.EQ. 'i')) THEN 
WRITE(11,*) ' ' 



ENDIF 

WRITE(2,*) NDP(LND(J,1),1), NDP(LND( J, 1) ,2) 
WRITE(2,*) ’ 999990 999990 ' 

DISSPLA PROGRAM GENERATION 

IF((CHAR3.EQ. ’D’).0R. (CHAR3.EQ. ’d')) 
WRITE(20,1010) NDP(LND(J,1),1), 
WRITE(20,1020) NDP(LND( J,2) , 1) , 
WRITE(20,1020) NDP(LND( J,3) , 1) , 
WRITE(20,1020) NDP(LND( J, 1) , 1) , 



THEN 

NDP(LND(J,1),2) 

NDP(LND(J,2),2) 

NDP(LND(J,3),2) 

NDP(LND(J,1),2) 



ENDIF 



CONTINUE 

FORMATC 1X,2(I3,2X),3X,3(F8. 5,2X)) 

F0RMAT(6X, 'CALL STRTPT( ' ,F8. 5, ' , ' ,F8. 5, ' )' ) 
FORMATC 6X, 'CALL CONNPT( ' ,F8. 5, ' , ' ,F8. 5, ' )' ) 



RETURN 

END 
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c 



c 

c 



40 



60 



SUBROUTINE ZERO 

ZERO A, B, C, AND P MATRICES 

COMPLEX A(50,50), B(50,50), C(50,50), P(50,100) 
INTEGER J, K, L 



C0MM0N/BLK8/A,B,C,P 

DO 50, J = 1, 50 

DO 40, K = 1, 50 

A(J,K) = CMPLX(0. 0,0. 0) 
B(J,K) = CMPLX(0. 0,0. 0) 
C( J,K) = CMPLX(0. 0,0. 0) 
CONTINUE 

DO 60, L = 1, 100 

P(J,L) = CMPLXCO. 0,0. 0) 
CONTINUE 



50 CONTINUE 
RETURN 
END 



SUBROUTINE BNDC ( I , BCOND , E 0 , SURBC , CHAR 1 , ALPHA , BETA , LINE , MODE , 
CX0RIGIN,Y0RIGIN,CHAR4) 

BOUNDARY CONDITION FILL FOR A PLANE WAVE 

EO = FIELD STRENGTH 

KO = WAVE NUMBER ( 2*P I /WAVELENGTH) 

BOUNDARY CONDITION FILL FOR A CYLINDRICAL BOUNDARY CONDITION 
EO = FIELD STRENGTH 

MODE = MODE NUMBER FOR CYLINDRICAL BOUNDARY CONDITIONS 

COMPLEX SURBC(IOO), VALUE(50), BCOND(IOO), ALPHA, BETA, LINE(50) 
COMPLEX K, RAK 

REAL NDP(200,2), EO, KR, KI , XORIGIN, YORIGIN, RA, RB 

INTEGER I, J, NDTOP, NDBOT, NDTOT, N0DES(200), PERND, BIND, MODE 

CHARACTER*! CHARI, CHAR4 

COMMON/BLK2/PERND,BIND 
COMMON/BLK3/NODES 
COMMON/ BLK5/NDP 

IF((CHAR1. EQ. 'P' ). OR. (CHARI. EQ. 'p' )) THEN 

IF((CHAR4. EQ. 'U').OR. (CHAR4.EQ. 'u')) THEN 
KR = REAL(CSQRT(BETA)) 

KI = ABS(AIMAG(CSQRT(BETA))) 

ELSE 

KR = 1.0 
KI = 0. 0 

ENDIF 

CHECK ON WHETHER WE ARE USING ER = ALPHA OR BETA 
IF(I.EQ. 1) THEN 

BCOND(l) = E0*EXP(KI*NDP(1,2))*CMPLX(C0S(KR*NDP(1,2)), 
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C SIN(KR*NDP(1,2))) 

BCOND(PERND) = EO*EXP(KI*NDP( 3,2) )*CMPLX(COS(KR*NDP( 3,2)), 

C SIN(KR*NDP(3,2))) 

BCOND(2) = EO*EXP(KI*NDP(7,2))*CMPLX(COS(KR*NDP(7,2)), 

C SIN(KR*NDP(7,2))) 

VALUE(l) = EO*EXP(KI*NDP(2,2))*CMPLX(COS(KR*NDP(2,2)) , 

C SIN(KR*NDP(2,2))) 

WRITE (10,*) BCOND(l) 

WRITE(10,1000) VALUE(l) 

WRITE(10,*) ' ' 

SURBC(l) = VALUE(l) 

ELSEIFd. LE. BIND) THEN 

NDTOP = NODES(I) + NODES( PERND-I+2) - 1 
NDBOT = NODES(I+l) + NODES(PERND-I+l) - 1 
NDTOT = NDTOP + NDBOT 

BCOND(PERND-I+l) = EO*EXP(KI*NDP(NDTOP+1 , 2) )*CMPLX(COS(KR* 

C NDP( NDTOP+1 , 2 ) ) , S I N( KR*NDP( NDTOP+ 1,2))) 

BCOND( I+l) = E0*EXP( KI*NDP( NDTOT, 2) )*CMPLX(COS(KR*NDP( NDTOT, 
C 2)) ,SIN(KR*NDP(NDTOT,2))) 

WRITE( 10,*) BCOND(PERND-I+2) 

DO 10, J = 2, NDTOP-1 

VALUE(J) = E0*EXP(KI*NDP(J,2))*CMPLX(C0S(KR*NDP(J,2)), 

C SIN(KR*NDP(J,2))) 

IF(J. EQ. 2) THEN 

SURBC( PERND-I+2) = VALUE(2) 

ELSEIF(J.EQ. (NDTOP-1)) THEN 

SURBC(I) = VALUE(NDTOP-l) 

ENDIF 

10 CONTINUE 

LINE(I) = VALUE(( NDTOP+1 )/2) 

WRITE(10,1000) (VALUE(J), J = 2, NDTOP-1) 

WRITECIO,*) BCOND(I) 

WRITEdO,*) ' ’ 

C 

ELSEIFd. EQ. BIND+1) THEN 

BC0ND(I+1) = EO*EXP(KI*NDP(6,2))*CMPLX(COS(KR*NDP(6,2)) , 

C SIN(KR*NDP(6,2))) 

VALUE(2) = E0*EXP(KI*NDP(2,2))*CMPLX(C0S(KR*NDP(2,2)), 

C SIN(KR*NDP(2,2))) 

VALUE(3) = E0*EXP(KI*NDP(3,2))*CMPLX(C0S(KR*NDP(3,2)), 

C SIN(KR*NDP(3,2))) 

VALUE(4) = E0*EXP(KI*NDP(4,2))*CMPLX(C0S(KR*NDP(4,2)), 

C SIN(KR*NDP(4,2))) 

WRITECIO,*) BC0ND(I+2) 

WRITECIO, 1000) VALUEC2), VALUE(3), VALUE(4) 

SURBCCBIND+3) = VALUE(2) 

SURBCC BIND+1) = VALUEC4) 

LINE(I) = VALUEC3) 

WRITECIO,*) BCOND(I) 

WRITECIO,*) ’ ' 

VALUEC2) = EO*EXP(KI*NDP(7,2))*CMPLX(COS(KR*NDP(7,2)) , 

C SIN(KR*NDP(7,2))) 

WRITECIO, 1000) VALUEC2) 

WRITECIO,*) BC0ND(I+1) 

SURBCC BIND+2) = VALUE(2) 
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END IF 
C 

ELSEIFC (CHARI. EQ. 'C').OR. (CHARI. EQ. 'c')) THEN 
C 

IF(I. EQ. 1) THEN 

RA = SQRT((NDP(2,1)-X0RIGIN)**2+(NDP(2,2)-Y0RIGIN)**2) 

RB = SQRT((NDP(1,1)-X0RIGIN)**2+(NDP(1,2)-Y0RIGIN)**2) 
WRITE(*,*) BETA 
K = CSQRT(BETA) 

RAK = RA*K 

CALL CYLBC( RA , RB , RAK , NDP( 1,1), NDP( 1,2), XORIGIN , YORIGIN, EO , 

C BCOND(l),SURBC(l),K) 

CALL CYLBC( RA , RB , RAK , NDP( 3 , 1 ) , NDP( 3 , 2 ) , XORIGIN , YORIGIN , EO , 

C BCOND(PERND),SURBC(PERND),K) 

CALL CYLBC( RA , RB , RAK , NDP( 7 , 1 ) , NDP( 7 , 2 ) , XORIGIN , YORIGIN , EO , 

C BCOND(2),SURBC(2),K) 

WRITEdO,*) BCOND(l), SURBC(l) 

WRITE(10,*) ' ' 

ELSEIF( I. LE. BIND) THEN 

NDTOP = NODES(I) + NODES ( PERND- 1+2) - 1 
NDBOT = NODES(I+l) + NODES(PERND-I+l) - 1 
NDTOT = NDTOP + NDBOT 

CALL CYLBC( RA , RB , RAK , NDP( NDTOP+1 , 1 ) , NDP( NDTOP+1 , 2 ) , XORIGIN , 
C YORIGIN, EO,BCOND(PERND-I+1),SURBC(PERND-I+1),K) 

CALL CYLBC( RA , RB , RAK , NDP( NDTOT , 1 ) , NDP( NDTOT , 2 ), XORIGIN , 

C YORIGIN,EO,BCOND(I+1) ,SURBC(I+1) ,K) 

WRITEdO,*) BCOND( PERND -1+2), SURBC( PERND -1+2) 

WRITEdO,*) BCOND(I), SURBC(I) 

WRITEdO,*) ' ' 

ELSEIFd. EQ. BIND+1) THEN 

WRITE(10,*) BCOND(I+2), SURBC(I+2) 

WRITEdO,*) BCOND(I), SURBC(I) 

WRITEdO,*) ’ ' 

CALL CYLBC( RA ,RB , RAK , NDP( 6 , 1 ) , NDP( 6 , 2 ) , XORIGIN , 

C YORIGIN, EO,BCOND( I+l) ,SURBC( I+l) ,K) 

WRITEdO,*) BCOND(I+l), SURBC(I+1) 

WRITEdO,*) ’ ' 

END IF 
C 

ELSE 

RETURN 

END IF 
C 

1000 F0RMAT(1X,50(E14. 8,2X,E14. 8,2X)) 

C 

RETURN 

END 



SUBROUTINE CYLBC ( RA , RB , RAK , X , Y , XORIGIN , YORIGIN , EO , BC , PSI , K) 

C 

COMPLEX SQRTMl, JORB, JIRB, JORAK, JIRAK, HORA, HIRA, HORB , HIRB 
COMPLEX DELTAN, AN, PSI, JORA, JIRA, K, RAK, BC 
REAL PI, XORIGIN, YORIGIN, X, Y, RA, RB , EO, PHI 
C 

PI = 4. 0*ATAN(1.0) 
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SQRTMl = CMPLXCO. 0,1. 0) 

PHI = ATAN2(X - XORIGIN, Y - YORIGIN) 

CALL BES1(CMPLX(RA,0. 0) , JORA, JIRA) 

CALL BES1(CMPLX(RB,0. 0) ,J0RB,J1RB) 

CALL BES1(RAK, J0RAK,J1RAK) 

CALL HAN1(CMPLX(RA,0. 0) ,H0RA,H1RA) 

CALL HAN1(CMPLX(RB,0. 0) ,H0RB,H1RB) 

C 

DELTAN = J1RB*(J1RAK*(H0RA-H1RA/RA)-K*(J0RAK-J1RAK/RAK)*H1RA)- 
CH1RB*(J1RAK*(J0RA-J1RA/RA)-K*(J0RAK-J1RAK/RAK)*J1RA) 

AN = -2. 0*SQRTM1/(PI*RA*DELTAN) 

BC = E0*COS(PHI) 

PSI = AN*J1RAK*BC 
C 

RETURN 

END 



SUBROUTINE BES1(Z, JO, Jl) 

Computing Bessel Functions for n = 0, 1 with 
Complex Argument Z. Direct Power Series Method for 
CABS(Z) . LE. 6 and Hankel's Asymptotic Formula for 
CABS(Z) .GT. 6. 

Written 11/5/87 by M. A. Morgan 
INTEGER M,M2 

REAL C(34) ,DM,F(34) ,G0,P(34) ,Pi,P2 

COMPLEX Z,Z2,Z3,Z4, J0,J1,AM,CL,P0,P1,Q0,Q1,C0,C1,S0,S1 
Pi=3. 1415927 
P2=2. 0/PI 

IF(CABS(Z). LE. 6. 0) THEN 

Utilizing the Direct Power Series Method 

G0= 1. 781072 
Z2=0. 5*Z 
CL=CLOG(GO*Z2) 

Computing F(m) = m ! and P(m) = 1 + 1/2 + 1/3 + ....+ 1/m 

F(l)=l. 0 
P(l)=l. 0 
DO 11 M=2,34 

F(M)=M*F(M-1) 

P(M)=P(M-1)+1. 0/M 
1 CONTINUE 

Computing Power Series Coefficients 

DM=-1. 0 
DO 22 M=l,34 

C(M)=DM/(F(M)*F(M)) 

DM=-DM 

22 CONTINUE 

C 
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C Computing JO and J1 

C 

J0=(1. ,0. ) 

J1=(0. ,0. ) 

M=0 

33 M=M+1 

M2=2*M 

AM=C(M)*(Z2**M2) 

J0=J0+AM 
J1=J1-M*AM 

IF((CABS(AM).GT. 1.0E-10).AND. (M. LT. 34)) GO TO 33 
J1=J1/Z2 
return 

ELSE 

Hankel' Asymptotic Formula (Abram. & Stegun p. 364) 

Z2=Z*Z 
Z3=Z*Z2 
Z4=Z-'Z3 

P0=1. 0-. 0703125/Z2+. 1121521/Z4 
Q0=-. 125/Z+. 0732422/Z3 
Pl=l. 0+. 1171875/Z2-. 1441956/Z4 
Ql=. 375/Z-. 10253906/Z3 
CO=CCOS(Z-. 25*PI) 

S0=CSIN(Z-. 25*PI) 

C1=CC0S(Z-. 75*PI) 

S1=CSIN(Z-. 75*PI) 

AM=CSQRT(P2/Z) 

J0=AM-'- ( P0*C0 -Q0*S0 ) 

Jl=AM--( P1*C 1 -Q1*S 1 ) 

END IF 
RETURN 
END 

SUBROUTINE HAN1( Z ,H0 ,H1) 



Computing Hankel Functions for n = 0, 1 with 
Complex Argument, Z. Direct Power Series Method for 
CABS(Z) . LE. 5 and Hankel's Asymptotic Formula for 
CABS(Z) .GT. 5. Written 11/6/87 by M. A. Morgan 

INTEGER M,M2 

REAL C(34),DM,F(34),G0,P(34),Pi,P2 

COMPLEX Z,Z2,Z3,Z4,J0,J1,Y0,Y1,AM,CL,P0,P1,Q0,Q1 

COMPLEX E0,E1,X0,X1,H0,H1, j 

PI=3. 1415927 

P2=2. 0/PI 

j=(0. .1. ) 

IF(CABS(Z). LE. 5. 0) THEN 

Direct Power Series Method 



G0= 1. 78072 
Z2=0. 5*Z 



no 



CL=CLOG(GO*Z2) 

C 

C Computing F(m) = m ! and P(m) = 1 + 1/2 + 1/3 + . . . . + 1/m 
C 

F(l)=l. 0 
P(l)=l. 0 
DO 11 M=2,34 

F(M)=M*F(M-1) 

P(M)=P(M-1)+1. 0/M 
11 CONTINUE 

C 

C Computing Power Series Coefficients 
C 

DM=-1. 0 
DO 22 M=l,34 

C(M)=DM/(F(M)*F(M)) 

DM=-DM 

22 CONTINUE 

C 

C Computing JO and J1 
C 

J0=(1. ,0. ) 

J1=(0. ,0. ) 

M=0 

33 M=M+1 

M2=2*M 

AM=C(M)*(Z2**M2) 

J0=J0+AM 

J1=J1-M*AM 

IF((CABS(AM).GT. 1.0E-10).AND. (M. LT. 34)) GO TO 33 
J1=J1/Z2 
C 

C Computing YO and Y1 
C 

M=0 

Y0=CL*J0 

Yl=Z2*CL*Jl-0. 5* JO 
44 M=M+1 

M2=2*M 

AM=C ( M ) *P ( M ) * ( Z 2**M2 ) 

Y0=Y0-AM 

Y1=Y1+M*AM 

IF((CABS(AM). GT. 1. OE-10). AND. (M. LT. 34)) GO TO 44 

Y0=P2*Y0 

Y1=P2*Y1/Z2 

H0=J0-j*Y0 

Hl=Jl-j*Yl 

RETURN 

ELSE 

C 

C Hankel' Asymptotic Formula (Abram. & Stegun p. 364 
C 

Z2=Z*Z 

Z3=Z*Z2 

Z4=Z*Z3 

P0=1. 0-. 0703125/Z2+. 1121521/Z4 
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Q0=-. 125/Z+. 0732422/Z3 
Pl=l. 0+. 1171875/Z2-. 1441956/Z4 
Ql=. 375/Z-. 10253906/Z3 
X0=(Z-. 25*PI) 

X1=(Z-. 75*PI) 

E0=CEXP(-j*X0) 

El=CEXP(-j*Xl) 

AM=CSQRT(P2/Z) 

H0=AM*(P0-j*Q0)*E0 

Hl=AM*(Pl-j*Qi)*Ei 

ENDIF 
RETURN 
END 

SUBROUTINE MARCH( I , IMX , NBMX , NA , NB , NC , INORM , SVEC ) 

THIS ROUTINE PERFORMS THE RICCATI TRANSFORM FIRST 
SWEEP, GENERATING AND STORING ON DISK 1 RMAT 
AND SVEC (FOR EACH MODE) AT EACH FORWARD STEP. 

COMPLEX RMAT(50,50),SVEC(50,100) 

COMPLEX A(50,50),B(50,50),C(50,50),P(50,100) 

COMPLEX D(50),SUM,DET 
REAL COND,DMAG 

INTEGER I , J , K , L , NA , NB , NC , PERND , B IND , IMX , INORM 
INTEGER NBMX 

CHARACTER*! CHAR2, CHAR3, CHAR4 

COMMON / BLK2/ PERND , B IND 
COMMON/BLK8/A,B,C,P 
COMMON/BLK9/CHAR2, CHAR3, CHAR4 

LOADING THEN INVERTING (B+A*RMAT) 

USING MINIMUM MEMORY SINGLE MATRIX TECHNIQUE 
DATE 1/29/80 FOR THIS CHANGE 

SKIPPING FIRST A*R (WHEN A = 0, ZERO R MATRIX) 

IF(I. EQ. 1) THEN 

DO 10, J = 1, NB 

DO 10, K = 1, NB 

RMAT(J,K) = (0. ,0. ) 

10 CONTINUE 
C RMAT = A*R 

ELSE 

C 

IF((CHAR2.EQ. ’l’).OR. (CHAR2.EQ. 'i')) THEN 
WRITE(19,*) ’ OLD R MATRIX’ 

DO 11, J = 1, 5 

WRITE(19,1000) (REAL(RMAT(J,K)), K = 1, 5) 

11 CONTINUE 
WRITE(19,*) ’ ’ 

WRITE(19,*) ' A MATRIX’ 

DO 12, J = 1, 5 

WRITE(19,1000) (REAL(A(J,K)), K = 1, 5) 

12 CONTINUE 
WRITE(19,*) ’ ’ 

ENDIF 

C 
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DO 30, K = 1, NB 

DO 20, J = 1, NB 

D(J) = (0. ,0. ) 

DO 20, L = 1, NA 

D(J) = D(J) + A(J,L)*RMAT(L,K) 

20 CONTINUE 

DO 30, J = 1, NB 

RMAT(J,K) = D(J) 

30 CONTINUE 

C 

IF((CHAR2. EQ. ' I ' ). OR. (CHAR2. EQ. 'i')) THEN 
WRITE(19,*) ' NEW R MATRIX’ 

DO 13, J = 1, 5 

WRITE( 19,1000) (REAL(RMAT(J,K)) , K= 1, 5) 

13 CONTINUE 
WRITE(19,*) ' ’ 

ENDIF 

C 

ENDIF 

C 

IF((CHAR2. EQ. ' l' ). OR. (CHAR2. EQ. ’ i’ )) THEN 
WRITE(19,*) ' B MATRIX' 

DO 14, J = 1, 5 

WRITE(19,1000) (REAL(B(J,K)) , K = 1, 5) 

14 CONTINUE 
WRITE(19,*) ' ' 

ENDIF 

C 

C RMAT = B + RMAT 

DO 40, J = 1, NB 

DO 40, K = 1, NB 

RMAT(J,K) = RMAT(J,K) + B(J,K) 

40 CONTINUE 

C 

IF((CHAR2.EQ. ’ I ’ ) . OR. ( CHAR2. EQ. 'i')) THEN 
WRITE(19,*) ' NEWEST R MATRIX' 

DO 15, J = 1, NB 

WRITE(9,1000) (REAL(RMAT(J,K)) , K = 1, NB) 
WRITE(19,1000) (REAL(RMAT(J,K)), K = 1, NB) 

15 CONTINUE 
WRITE (9,*) ’ ' 

WRITE(19,*) ' ' 

ENDIF 

C INVERTING THE MATRIX (B + A*R) 

CALL CSMINVC RMAT , NBMX , NB , DET , COND , INORM) 

C 

DMAG = CABS(DET) 

WRITE(*,*) I, NBMX, NB, DMAG, COND 
C 

C COMPUTING THE NEW S -VECTORS 
C 

C SKIPPING FIRST A*S (A = 0) 

IF (I.EQ. 1) THEN 
CONTINUE 

ELSE 

DO 70, K = 1, PERND 
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DO 60, J = 1, NB 

D(J) = (0. 0,0. 0) 

DO 60, L = 1, NA 

D(J) = D(J) + A(J,L)*SVEC(L,K) 

60 CONTINUE 

DO 70, J = 1, NB 

P(J,K) = P(J,K) - D(J) 

70 CONTINUE 

ENDIF 

C FINAL SVEC MULTIPLICATION 

DO 100, K = 1, PERND 
DO 90, J = 1, NB 

D( J) = (0. 0,0. 0) 

DO 90, L = 1, NB 

D(J) = D(J) + RMAT(J,L)*P(L,K) 

90 CONTINUE 

DO 100, J = 1, NB 

SVEC(J,K) = D(J) 

100 CONTINUE 

C 

C STORING I+l SVEC ON DISK 7 

DO 110 J = 1, NB 

WRITE(7) (SVEC(J,K), K = 1, PERND) 

110 CONTINUE 
C 

C FINAL RMAT MULTIPLICATION 

IFCI.EQ. IMX) RETURN 
DO 130, J = 1, NB 

DO 120, K = 1, NC 

D(K) = (0. 0,0.0) 

DO 120, L = 1, NB 

D(K) = D(K) - RMAT(J,L)*C(L,K) 

120 CONTINUE 

DO 130, K = 1, NC 

RMAT(J,K) = D(K) 

130 CONTINUE 
C 

C STORING I+l RMAT ON DISK 7 
DO 140, J = 1, NB 
WRITEC7) (RMAT(J,K), K = 1, NC) 

140 CONTINUE 
C 

1000 F0RMAT(10(E9. 3,1X)) 

C 

RETURN 

END 

C 

C 

SUBROUTINE SWEEP( IMX , NABC , SURBC , LINE , CHARI , U , BCOND , PS I , ANS , CHARS ) 
C THIS ROUTINE PERFORMS THE RICCATI TRANSFORM BACKSWEEP 

C FROM I=IMX TO 1=1, RECALLING RMAT AND 

C SVEC FROM DISK 7 AT EACH BACKSTEP TO FORM 

C THE NODE VECTORS PSI, THEN STORING THESE ON 

C DISK 8 FOR EACH APPLIED DIRICHLET B. C. 

COMPLEX RMAT(50,100) ,PSI(50 , 100) ,SVEC( 50 , 100) ,D( 100) ,TEMP 
COMPLEX SURBC( 100) ,ANS( 100) ,LINE(50) ,ANSB(50) ,U( 100,100) 
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c 

c 



c 

c 



10 



15 

20 

C 

C 



30 

C 



40 

C 



50 



60 

C 



70 



COMPLEX BCOND(IOO) 

REAL ERROR? , ERRORD , TERRN , TERRD , ATSERR 
INTEGER I,J,K,L,NABC(100,3) ,IMX,PERND,BIND 
CHARACTER*! CHARI, CHAR5 

C0MM0N/BLK2/PERND,BIND 

IF((CHAR5.EQ. 'M' ). OR. (CHAR5. EQ. 'm')) THEN 
RETURN 

END IF 



INITIAL DISK READ AT IMX (R = 0, NOT WRITTEN, => S IS READ FIRST) 
WRITE (*,*) ' ' 

DO 90, I = IMX, 1, -1 

WRITE(*,1050) (IMX-I+1), IMX 
IF(I.EQ. IMX) THEN 

DO 10, J = NABC(I,2) ,1,-1 

RAPlir^PArF 7 

READ(7) (PSI(J,K), K = 1, PERND) 

BACKSPACE 7 
CONTINUE 

DO 20, J = 1, NABC(I,2) 

DO 15, K = 1, PERND 

U(IMX,K) = PSI(J,K) 

CONTINUE 

CONTINUE 

SUBSEQUENT DISK READS 
ELSE 

READ R MATRIX 

DO 30, J = NABC(I,2), 1, -1 
BACKSPACE 7 

READ(7) (RMAT(J,K), K = 1, NABC(I,3)) 

BACKSPACE 7 
CONTINUE 
READ S VECTOR 

DO 40, J = NABC(I,2), 1, -1 
BACKSPACE 7 

READ(7) (SVEC(J,K), K = 1, PERND) 

BACKSPACE 7 
CONTINUE 

MULTIPLY RMAT = RMAT*PSI 
DO 60, J = 1, NABC(I,2) 

DO 50, K = 1, PERND 
D(K) = (0. 0,0. 0) 

DO 50, L = 1, NABC(I,3) 

D(K) = D(K) + RMAT(J,L)*PSI(L,K) 

CONTINUE 

DO 60, K = 1, PERND 
RMAT(J,K) = D(K) 

CONTINUE 

PSI = RMAT + SVEC 
DO 70, J = 1, NABC(I,2) 

DO 70, K = 1, PERND 

PSI(J,K) = RMAT(J,K) + SVEC(J,K) 

CONTINUE 

ANSB(I) = (0. 0,0.0) 
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75 

80 

90 

C 

94 

98 

C 



100 



C 



110 



C 



DO 80, J = 1, NABC(I,2) 

DO 75, K = 1, PERND 

IF((J. EQ. NABC(I,2)).0R. (I.EQ. 1)) THEN 
U(I,K) = PSI(J,K) 

ELSEIFC J. EQ. 1) THEN 

U(2*IMX-I,K) = PSI(J,K) 

ELSEIFC J. EQ. (NABC(I,2)+l)/2) THEN 

ANSB(I) = ANSB(I) + PSI( J,K)*BCOND(K) 

END IF 
CONTINUE 
CONTINUE 

ENDIF 
CONTINUE 
WRITEC8,*) ' ' 

DO 98, J = 1, PERND 

ANS(J) = (0. 0,0. 0) 

DO 94, L = 1, PERND 

ANS(J) = ANS(J) + U(J,L)*BCOND(L) 

CONTINUE 

CONTINUE 

TERRN = 0. 0 
TERRD = 0. 0 
DO 100, I = 1, PERND 

ERRORP = (CABS(ANSCI) - SURBC( I ) ) )**2 
ERRORD = (CABS(SURBC(I)))*-''2 
TERRN = TERRN + ERRORP 
TERRD = TERRD + ERRORD 

WRITE(13,1020) I, BCOND(I), SURBC(I), ANS(I), ERRORP 
CONTINUE 
WRITEC13,*) ' ' 

ATSERR = (TERRN/TERRD)**0. 5 
WRITE(13,1030) ATSERR 

WRITE(-'-,*) 'RMS ERROR (FOR THE PERIMETER) = ', ATSERR 
WRITEC13,*) ' ’ 

IFCCCHARl. EQ. ’C').0R. (CHARI. EQ. 'c')) THEN 
RETURN 

ELSE 

TERRN = 0. 0 
TERRD = 0.0 
DO 110, 1=2, BIND+1 

ERRORP = (CABS(ANSBd) - LINE(I)))**2 
ERRORD = (CABS(LINE(I)))**2 
TERRN = TERRN + ERRORP 
TERRD = TERRD + ERRORD 

WRITE(13,1025) (I-l), LINE(I), ANSB(I), ERRORP 
CONTINUE 
WRITE(13,*) ’ ' 

ATSERR = (TERRN/TERRD)**0. 5 
WRITE(13,1040) ATSERR 

WRITE(*,*) 'RMS ERRROR (FOR BISECTION SEGMENT) = ', ATSERR 

ENDIF 

DO 120, 1=1, PERND 
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WRITE(30, 1060) (U(I,J), J = 1, PERND) 

120 CONTINUE 
C 

1020 F0RMAT(1X,I3,2X,3(E14. 8,1X,E14. 8,3X) ,F10. 6) 

1025 F0RMAT(1X,I3,2X,4(E14. 8,2X) ,F10. 6) 

1030 FORMATdX,’ RMS ERROR (FOR THE PERIMETER) = ' ,F12. 6) 

1040 FORMATdX,' RMS ERROR (FOR BISECTION SEGMENT) = ' ,F12. 6) 

1050 FORMATdX, 13,' TOTAL BACKSWEEP ROWS OUT OF ',13,' COMPLETED') 
1060 FORMATdX, 100(E8. 2, E8. 2,2X)) 

C 

RETURN 

END 

C 

CC 

C 

SUBROUTINE CSMINV( A , NDIM , N , DETERM , COND , INORM) 

C 

C INORM - FLAG TO NORMALIZE COLUMNS AND ROWS OF MATRIX A 
C MATRIX NORMALIZATION BY M. A. MORGAN 

C APRIL 24,1978 

C 

C A 

C NDIM 

C N 

C DETERM 

C COND 

C INORM 

C 
C 

COMPLE XA(50, 50), PI V0T(50), AMAX, T, SWAP, DETERM, U 

INTEGER I , J , K , L, IPI VOT( 50) , INDEX( 50 , 2) , IROW , ICOLUM , LI , JROW 

INTEGER JCOLUM,N, INORM 

REAL TEMP , ALPHA( 50 ) , COL( 50 ) , ROW( 50) , AJK , SUMAXA , SUMROW , SUMAXI 
C 

IF(NDIM. GT. 50) THEN 

WRITE(*,*) ' ERROR IN INVERTION CALL. . . DIMENSION > 50 ’ 
STOP 

END IF 

IF(N. GT. NDIM) THEN 

WRITE(*,*) ' ERROR IN INVERTION CALL. . . N > MAX DIM. ' 
STOP 

END IF 

IF( INORM. NE. 1) GO TO 7 
DO 3 K = I, N 

COL(K) =0.0 
DO 1 J = I,N 

AJK = CABS(A(J,K)) 

IF(AJK. GT. COL(K)) COL(K) = AJK 

1 CONTINUE 

DO 2 J = 1, N 

2 A(J,K) = A( J,K)/COL(K) 

3 CONTINUE 

C ROW NORMALIZING 

DO 6 J = 1, N 

ROW(J) =0.0 
DO 4 K = I, N 



- MATRIX TO INVERT 



DETERMINATE OF A 
CONDITION NUMBER OF A 
INTEGER NORMALIZATION FLAG 



INPUT/OUTPUT 

INPUT 

INPUT 

OUTPUT 

OUTPUT 

INPUT 
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AJK = CABS(A(J,K)) 

IF(AJK. GT. ROW(J)) ROW(J) = AJK 

4 CONTINUE 

DO 5 K = 1, N 

5 A(J,K) = A(J,K)/ROV(J) 

6 CONTINUE 

7 CONTINUE 

DETERM = CMPLX(1. 0,0. 0) 

SUMAXA = 0. 0 
DO 20 J = 1, N 

ALPHA(J) =0.0 
SUMROW = 0.0 
DO 10 I = 1, N 

ALPHA(J) = ALPHA(J) + A(J,I)* CONJG( A( J, I ) ) 

10 SUMROW = SUMROW + CABS(A(J,I)) 

ALPHA(J) = SQRT(ALPHA(J)) 

IF( SUMROW. GT. SUMAXA) SUMAXA = SUMROW 
20 IPIVOT(J) = 0 

DO 600 I = 1, N 
C 
C 

AMAX = CMPLX(0. 0,0. 0) 

DO 105 J = 1, N 

IF (IPIVOT(J)-l) 60, 105, 60 
60 DO 100 K = 1, N 

IF (IPIVOT(K)-l) 80, 100, 740 

80 TEMP = AMAX*CONJG(AMAX) - A( J,K)*CONJG( A( J,K) ) 

IF(TEMP)85,85,100 

85 IROW = J 

ICOLUM = K 
AMAX = A(J,K) 

100 CONTINUE 

105 CONTINUE 

IPIVOT(ICOLUM) = IPIVOT(ICOLUM) + 1 
C 
C 

IF (IROW- ICOLUM) 140, 260, 140 
140 DETERM = -DETERM 

DO 200 L = 1, N 

SWAP = A(IROW,L) 

A(IROW,L) = A( ICOLUM, L) 

200 A( ICOLUM, L) = SWAP 

SWAP = ALPHA(IROW) 

ALPHA(IROW) = ALPHA( ICOLUM) 

ALPHA( ICOLUM) = SWAP 
260 INDEX(I,1) = IROW 

INDEX(I,2) = ICOLUM 
PIVOT(I) = A( ICOLUM, ICOLUM) 

U = PIVOT(I) 

DETERM = DETERM--U 

DETERM = DETERM/ALPHAC ICOLUM) 

TEMP = PIVOT(I)*CONJG(PIVOT(I)) 

IF(TEMP)330, 720,330 
C 
C 

330 A( ICOLUM, ICOLUM) = CMPLX( 1. 0 ,0. 0) 
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350 

C 

C 

380 

400 

450 

550 

600 

C 

C 

620 

630 

705 

710 

900 

910 

950 

955 

720 

730 

740 

C 

C 

C 

C 

C 



DO 350 L = 1, N 
U = PIVOT(I) 

A(ICOLUM,L) = A(ICOLUM,L)/U 



DO 550 LI = 1, N 

IF(Ll-ICOLUM) 400, 550, 400 
T = A(L1,IC0LUM) 

A(L1,IC0LUM) = CMPLX(0. 0,0. 0) 
DO 450 L = 1, N 

U = A(ICOLUM,L) 

A(L1,L) = A(L1,L) - U*T 

CONTINUE 

CONTINUE 



DO 710 I = 1, N 

L = N + 1 - I 

IF (INDEX(L,1) - INDEX(L,2)) 630, 710, 630 
JROW = INDEX(L,1) 

JCOLUM = INDEX(L,2) 

DO 705 K = 1, N 

SWAP = A(K,JROW) 

A(K,JROW) = A(K, JCOLUM) 

A(K, JCOLUM) = SWAP 
CO.NTINUE 

CONTINUE 
SUMAXI = 0. 0 
DO 910 I = 1, N 

SUMROW =0.0 

DO 900 J = 1, N 

SUMROW = SUMROW + CABS(A(I,J)) 

IF( SUMROW. GT. SUMAXI) SUMAXI = SUMROW 
CONTINUE 

COND = SUMAXA*SUMAXI 
IFdNORM. NE. 1) GO TO 955 
DO 950 K = 1, N 

DO 950 J = 1, N 

A(J,K) = A(J,K)/(ROW(K)*COL(J)) 

CONTINUE 
RETURN 
WRITE (*,730) 

FORMATC ' MATRIX IS SINGULAR’) 

RETURN 

END 

SUBROUTINE S AVE ( BCOND , ANS , U , OFFSET , PERND , CHAR5 , DPER , K , XORG , YORG , 
CNRE S , MRE S , LB I AS , GB I AS , MAXD ) 

THIS SUBROUTINE SAVES THE ESSENCE OF THE FINITE ELEMENT PROBLEM 
TO A DATA FILE CALLED " F3. DAT ". tHIS DATA IS NECESSARY TO 
SOLVE THE FIELD FEEDBACK FORMULATION, (F3). 

COMPLEX BCOND( 100) ,ANS( 100) ,U( 100 , 100) 

RE AL OFFSET , MESH( 0 : 200 , 5 ) , DPER , K , XORG , YORG , MRE S , MAXD 
INTEGER PERND , I , J , NRES , LBI AS , GBIAS 
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CHARACTER*! CHARS 
C 

COMMON/BLKl/MESH 

C 

IF((CHAR5. EQ. 'M' ). OR. (CHARS. EQ. 'm' )) THEN 
RETURN 

ENDIF 

C 

WRITE(40,*) PERND 
WRITE(40,*) OFFSET 
WRITE(40,*) DPER 
WRITE (40,*) K 
WRITE (40,*) XORG 
WRITE (40,*) YORG 
WRITE(40,*) NRES 
WRITE(40,*) MRES 
WRITE(40,*) LBIAS 
WRITE(40,*) GBIAS 
WRITE (40,*) MAXD 
C 

DO 10, I = 1, 4 

DO 10, J = 1, PERND 
WRITE(40,*) MESH(J,I) 

10 CONTINUE 
C 

DO 20, J = 1, PERND 

WRITE(40,*) BCOND(J) 

20 CONTINUE 

C 

DO 30, J = 1, PERND 

WRITE(40,*) ANS(J) 

30 CONTINUE 

C 

DO 40, I = 1, PERND 

DO 40, J = 1, PERND 

WRITE(40,*) U(I,J) 

40 CONTINUE 

C 

RETURN 

END 
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APPENDIX D. VARINT CONVERGENCE PROGRAM 



TEST OF VARINT CONVERGENCE 

COMPLEX F(3,3), ALPHA, BETA, EXACT, CENTER, LEFT, RIGHT, TOP 

COMPLEX BOTTOM, SUM, CALC 

REAL X(3), Y(3), D, AREA, KR, PI, ERROR 

INTEGER I 

OPEN (UNIT = 1, FILE = 'C: MSFORT TEST.DAT’, STATUS = 'UNKNOWN') 
ALPHA = (1. 0,0. 0) 

BETA = (1. 0,0. 0) 

PI = 4. 0*ATAN(1. 0) 

DO 5, 1=1, 100 

D = 2. 0*PI-'-FL0AT(I)/100. 0 
KR = REAL(CSQRT(BETA)) 

X(l) = 0. 0 
Y(l) = 0.0 
X(2) = D 
Y(2) = 0.0 
X(3) = 0. 0 
Y(3) = D 

CALL VARINTC X , Y , F , ALPHA , BETA , ARE A ) 

EXACT = CMPLX(C0S(KR*0. 0) ,SIN(KR*0. 0)) 

CENTER = 4. 0*F(1,1) 

RIGHT = -2.0*F(l,2)*CMPLX(COS(KR*Y(2)),SIN(KR*Y(2))) 

TOP = -2. 0*F(1,2)*CMPLX(C0S(KR*Y(3)),SIN(KR*Y(3))) 

LEFT = -2. 0*F( l,2)*CMPLX(COS(KR*Y(2)) ,SIN(KR*Y(2))) 

BOTTOM = -2. 0*F(1,2)*CMPLX(COS(-KR*Y(3)),SIN(-KR*Y(3))) 

SUM = TOP + BOTTOM + LEFT + RIGHT 
CALC = SUM/ CENTER 
ERROR = (CALC -EXACT) /EXACT 
WRITE( 1, 1000) I ,D, EXACT, CALC, ERROR 
CONTINUE 

CLOSE(l) 

F0RMAT( IX, 13, 1X,F8. 5, 1X,2(F8. 5, 1X,F8. 5, IX) ,E13. 6) 

STOP 

END 



SUBROUTINE VARINT( X , Y , F , ALPHA , BETA , AREA) 

GENERATING VARIATIONAL FINITE ELEMENT AREA INTEGRATIONS OF THE 
LINEAR BASIS FUNCTION LAGRANGIAN FOR THE HELMHOLTZ EQUATION. 
THESE ARE RETURNED IN F(3,3). X(3) AND Y(3) ARE THE WAVENUMBER 
NORMALIZED CARTESIAN COORDINATES OF THE TRIANGLE VERTICES. 

X = Ko*x, Y = Ko*y ALPHA AND BETA ARE COMPLEX 
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MATERIAL PARAMETERS WITHIN THE ELEMENT. 

FOR TM INCIDENCE: ALPHA = 1/ur; BETA = er 
FOR TE INCIDENCE: ALPHA = 1/er; BETA = ur 

OUTPUTS : F(3,3) - FINITE ELEMENT AREA INTEGRATION 

AREA - AREA OF A TRIANGLE 

COMPLEX ALPHA, BETA, B12, F(3,3) 

REAL X(3), Y(3), T(3,3), AREA, DET 
INTEGER L, K 

DET = ABS(X(2)*Y(3) + X(3)*Y(1) + X(1)*Y(2) - X(3)*Y(2) - 
CX(1)*Y(3) - X(2)*Y(1)) 

AREA = ABS(0. 5*DET) 

B12 = BETA/12. 

T(l,l) = (Y(2) - Y(3))/DET 

T(l,2) = (Y(3) - Y(1))/DET 

T(l,3) = (Y(l) - Y(2))/DET 

T(2,l) = (X(3) - X(2))/DET 

T(2,2) = (X(l) - X(3))/DET 

T(2,3) = (X(2) - X(1))/DET 

T(3,l) = (X(2)*Y(3) - X(3)*Y(2))/DET 

T(3,2) = (X(3)*Y(1) - X(1)*Y(3))/DET 

T(3,3) = (X(1)*Y(2) - X(2)*Y(1))/DET 

DO 10, K = 1, 3 

DO 10 L = 1, 3 

F(K,L) = ALPHA*(T(1,K)*T(1,L) + T( 2 ,K)*T( 2 , L) ) - B12 
IF(K. EQ. L) F(K,L) = F(K,L) - B12 
F(K,L) = AREA*F(K,L) 

RETURN 

END 
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APPENDIX E. FIELD FEEDBACK PROGRAM 

c 

C FIELD FEEDBACK FORMULATION PROGRAM 

C WRITTEN BY T. B. WELCH 

C 

C w/ PROGRAMMING IDEAS FROM PROF M. A. MORGAN 
C 

C BCOND 

C OFFSET 

C ANS 

C DPER 

C 

C LB I AS 

C GBIAS 

C MAXD 

C 

C K 

C XORG 

C YORG 

C NRES 

C 

C U 

C 
C 

C PERND 

C MESH 

C 
C 
C 
C 
C 
C 
C 
C 

C T 

C 

C CNVEC 

C MRES 

C 
C 

COMPLEX BCOND( 100) ,ANS( 100) ,U( 100 , 100) ,T( 100 , 100) ,CNVEC( 100) 

REAL OFFSET , ME SH ( 1 0 0 , 8 ) , DPER , K , XORG , YORG , MRE S , MAXD 
INTEGER PERND , I , J , NRE S , LB I AS , GB I AS 
C 

OPEN (UNIT = 40, FILE = 'C: MSFORT F 3. DAT' , STATU S=' UNKNOWN' ) 

OPEN (UNIT = 50, FILE = 'C: MSFORT FFPAT. DAT' ,STATUS=' UNKNOWN' ) 

C 

CALL INPUT( BCOND , ANS , U , OFFSET , PERND , ME SH , DPER , K , NRES , XORG , YORG , 
CMRES,LBIAS, GBIAS, MAXD) 

CALL TMAT( U , PERND , ME SH , T , OFFSET , DPER , BCOND , MRE S , LB I AS , GB I AS , MAXD ) 
CALL CNSOLV(T, BCOND, CNVEC, PERND) 



- BOUNDARY CONDITIONS 

- OFFSET IN WAVELENGTHS (PERIMETER TO BOUNDARY) 

- CALCULATED PS I VALUES ON PERIMETER 

- DESIRED PERCENT ERROR SCALE FACTOR FOR GREEN'S 
FUNCTION INTEGRAL PATCH STEPPING 

- BIAS THAT IS ADDED TO THE GFI STEP FOR < 1. 0 

- BIAS THAT IS ADDED TO THE GFI STEP FOR >1.0 

- MAXIMUM DISTANCE BEYOND WHICH NO CONTRIBUTION IS MADE 
TO THE GFI 

- WAVENUMBER 

- X ORIGIN 

- Y ORIGIN 

- NUMBER OF EVENLY SPACED POINTS DESIRED FOR THE FAR 
FIELD CALCULATIONS (360/NRES = ANGULAR RESOLUTION) 

- MATRIX THAT RELATES EACH BOUNDARY NODE VALUE TO THE 
UNKNOWN PERIMETER NODE VALUE. MULTIPLY U BY A DRIVING 
VECTOR (ON BOUNDARY) TO FIND PERIMETER VALUES. 

- NUMBER OF PERIMETER NODES 

- GEOMETRY ARRAY CONTAINING: 

1 - X POSITION OF PERIMETER NODES 

2 - Y POSITION OF PERIMETER NODES 

3 - X UNIT NORMAL OF PERIMETER NODES 

4 - Y UNIT NORMAL OF PERIMETER NODES 

5 - X POSITION OF GEOMETRIC CONTOUR NODES 

6 - Y POSITION OF GEOMETRIC CONTOUR NODES 

7 - X POSITION OF BOUNDARY NODES 

8 - Y POSITION OF BOUNDARY NODES 

- MATRIX THAT RELATES PERIMETER VALUES BACK OUT 
TO THE BOUNDARY VIA A GREEN'S FUNCTION INTEGRAL 

- VECTOR OF SCATTERED FIELD BACK ONTO THE BOUNDARY 

- MESH RESOLUTION 
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CALL FFLD ( CNVEC , PERND , ME SH , U , OFF SET , K , NRE S , XORG , YORG ) 
C 

CLOSE (40) 

CLOSE(SO) 

C 

STOP 

END 



SUBROUTINE INPUT( BCOND , ANS , U , OFFSET , PERND , MESH , DPER , K , NRES , XORG , 
CYORG , MRES , LBI AS , GB I AS , MAXD) 

THIS SUBROUTINE READS THE FINITE ELEMENT PROBLEM DATA FROM 
THE DATA FILE CALLED " F3. DAT THIS DATA IS NECESSARY TO 
SOLVE THE FIELD FEEDBACK FORMULATION, (F3). 

COMPLEX BCOND( 100) ,ANS( 100) ,U( 100 , 100) 

REAL OFFSET, MESH( 100, 8), DPER, K, XORG, YORG, MRES, MAXD 
INTEGER PERND , I , J ,NRES , LBIAS , GBIAS 
C 

WRITE(*,Vc) ' reading INPUT DATA ' 

C 

RE AD (40,*) PERND 
READ(40,*) OFFSET 
READ(40,*) DPER 
READ(40,*) K 
RE AD (40,*) XORG 
READ(40,*) YORG 
READ(40,*) NRES 
READ(40,*) MRES 
READ(40,*) LBIAS 
READ(40,*) GBIAS 
READ(40,*) MAXD 
C 

WRITE(6,*) ’ NUMBER OF PERIMETER NODES 
WRITE (6,*) ’ BOUNDARY CONTOUR OFFSET 
WRITE(6,*) ' DESIRED GFI SCALE FACTOR 
WRITE(6,*) ' GFI STEP BIAS FOR <1.0 

WRITE(6,*) ' GFI STEP BIAS FOR >1.0 

WRITE(6,*) ' MAX DIST > NO CONTRIBUTION 
WRITE(6,*) ' WAVENUMBER 
WRITE(6,*) ' X ORIGIN 

WRITE(6,*) ' Y ORIGIN 

WRITE(6,*) ' NUMBER OF NODES FOR SIGMA 
WRITE(6,*) ' REQUESTED MESH RESOLUTION 
C 

DO 10, I = 1, 4 

DO 10, J = 1, PERND 

READ(40,*) MESH(J,I) 

10 CONTINUE 

C 

DO 20, J = 1, PERND 

READ(40,*) BCOND(J) 

20 CONTINUE 

C 

DO 30, J = 1, PERND 



= ' , PERND 
= ' , OFFSET 
= ' ,DPER 
= ' , LBIAS 
= ' , GBIAS 
GFI = ' ,MAXD 
= ',K 
= ' ,XORG 
= ' ,YORG 
= ' ,NRES 
= ' ,MRES 
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30 

C 



40 

C 



50 

C 



READ(40,*) ANS(J) 
CONTINUE 

DO 40, 1=1, PERND 

DO 40, J = 1, PERND 

READ(40,*) U(I,J) 

CONTINUE 



DO 50, 1=1, PERND 

MESH(I,5) = MESH(I,1) 
MESH(I,6) = MESH(I,2) 
MESH(I,7) = MESH(I,1) 
MESH(I,8) = MESH(I,2) 
CONTINUE 



+ MESH(I,3)*OFFSET/2. 0 
+ MESH(I,4)*0FFSET/2. 0 
+ MESH(I,3)*OFFSET 
+ MESH(I,4)*0FFSET 



RETURN 

END 



SUBROUTINE TMT( U , PERND , MESH , T , OFFSET , DPER , BCOND , MRES , LBI AS , GBIAS , 
CMAXD) 

THIS SUBROUTINE CALCULATES THE GREEN’S FUNCTION INTEGRAL 
(FOR A SINGLE BASIS FUNCTION BOUNDARY CONDITION) GEOMETRIC 
PERIMETER WITH RESPECT TO EACH OF THE OFFSET BOUNDARY NODES. 

THE INTEGRATION IS REPEATED UNTIL EACH BOUNDARY CONDITION HAS BEEN 
INDIVIDUALLY APPLIED AND INTEGRATED WITH RESPECT TO EACH OFFSET 
BOUNDARY NODE. THESE VALUES ARE RETURNED IN THE " T ” MATRIX 
FOR USE IN THE FIELD FEEDBACK FORMULATION. THE MATRIX IS 
ORGANIZED, T m,n. 

COMPLEX U( 100,100) ,T( 100, 100) ,PVEC( 100) ,PDVEC( 100) 

COMPLEX J,H0RP1 ,H0RP2,H1RP1 ,SUM,TEMP( 100) ,BCOND( 100) 

COMPLEX H1RP2 , PS I , PS IRP , INTEGRAL , DPS IC 
REAL RMRP , NORM 1 1 , NORM 12 , DOT , DPER , D I STM , MAXD 
REAL OFFSET, MESH( 100,8) ,DIST,R,DZ,DL,MRES 

INTEGER I , M , N , NN , STEP , FN , SN , PERND , MM , STEPMX , STEPMN , LB I AS , GB I AS 
OPEN (UNIT = 2, FILE = 'C: MSFORT TMAT. DAT’ , STATUS = ’UNKNOWN’) 

J = (0. 0,1.0) 

STEPMX = INT(DPER*(-48. 0*0FFSET + 17.2) + LBIAS) 

STEPMN = INT(DPER + GBIAS) 

WRITE(*,*) ’ LOADING T MATRIX ’ 

WRITE(*,*) ’ MAXIMUM STEP = ’, STEPMX,’, MINIMUM STEP = ’, STEPMN 
WRITE(*,*) ’ MAXIMUM DISTANCE FOR ANY CONTRIBUTION = ’ ,MAXD 
DO 40, M = 1, PERND 

DO 5, 1=1, PERND 

IF(M. EQ. I) THEN 

PVEC(I) = (1.0 + U(I,M))/2.0 
PDVEC(I) = (1.0 - U( I, M)) /OFFSET 

ELSE 

PVEC(I) = U(I,M)/2.0 
PDVEC(I) = -U( I, M) /OFFSET 
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5 

C 



C 

C 



C 



C 



ENDIF 

CONTINUE 



DO 30, N = 1, PERND 

WRITE(*,1000) M, N, PERND 
SUM = (0. 0,0. 0) 

DO 20, NN = 1, PERND 
FN = NN 

IF(NN. EQ. PERND) THEN 
SN = 1 

ELSE 

SN = NN + 1 

ENDIF 

DIST = SQRT((MESH(FN,5)-MESH(SN,5))**2+(MESH(FN,6) 

C -MESH(SN,6))**2) 

INTEGRAL = (0. 0,0. 0) 

DISTM = SQRT((MESH(N,7)-MESH(FN,5))**2+(MESH(N,8)-MESH(FN,6))**2) 

R = MESH(SN,5) - MESH(FN,5) 

DZ = MESH(SN,6) - MESH(FN,6) 

DL = SQRT(R**2 + DZ**2) 

NORMll = -DZ/DL 
NORM12 = R/DL 

IF(DISTM. GT. MAXD) THEN 
GOTO 20 

ELSEIFC DISTM. LE. 1. 0) THEN 
STEP = STEPMX 

ELSE 

STEP = STEPMN 

ENDIF 

IF(STEP. LT. 1) STEP = 1 

DO 10, 1=1, STEP+1 
IFCI.EQ. 1) THEN 

RMRP = SQRT((MESH(N,7)-(MESH(FN,5)+0. 25*(MESH(SN,5)- 
MESH( FN , 5 ) ) /FLOAT( STEP ) ) )**2+( MESH( N , 8 ) - ( MESH( FN , 6 )+ 

0. 25*( MESH( SN , 6) -MESH( FN , 6) ) /FLOAT( STEP) ) )**2 ) 

DPSIC = PDVEC(FN) + 0. 25*(PDVEC(SN) -PDVEC(FN) )/STEP 
PSIRP = PVEC(FN) + 0. 25*(PVEC(SN)-PVEC(FN))/STEP 
DOT = (N0RM11*(MESH(N,7) - (MESH(FN,5)+0.25*(MESH(SN,5)- 
MESH(FN,5))/STEP))+(NORM12*(MESH(N,8)-(MESH(FN,6)+0. 25* 
(MESH(SN,6)-MESH(FN,6))/STEP))))/RMRP 
ELSEIFC I. EQ. STEP+1) THEN 

RMRP = SQRT((MESH(N,7)-(MESH(FN,5)+(FL0AT(STEP)-0. 25)*( 
MESHC SN , 5 ) -MESHC FN , 5 ) ) /FLOATC STEP) ) )**2+( MESH( N , 8 ) - ( MESH 
( FN , 6 ) +( FLOATC STEP ) - 0 . 25 ) *C MESHC SN , 6 ) -MESHC FN , 6 ) ) /FLOAT 
CSTEP)))**2) 

DPS I C=PDVEC C FN ) +C FLOATC STEP)-0.25)*C PDVEC C SN ) - PDVEC C FN ) ) 
/CSTEP) 

PSIRP = PVECCFN)+CFLOATCSTEP)-0. 25)*CPVECCSN)-PVECCFN))/ 
CSTEP) 

DOT = CNORM11*CMESHCN,7)-CMESHCFN,5)+CFLOATCSTEP)-0.25)* 
C MESHC SN , 5 ) -MESHC FN , 5 ))/ STEP) ) +C NORM12*C MESHC N , 8 ) - C MESH 
CFN, 6)+C FLOATC STEP) -0. 25)*CMESHC SN, 6) -MESHCFN, 6) )/STEP) 
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c 

10 

20 

30 

40 

C 



50 

55 

C 



60 

C 

1000 

C 

C 



C 

C 

C 

C 

C 

C 

C 

C 



C )))/RMRP 

ELSE 

RMRP = SQRT((MESH(N,7)-(MESH(FN,5)+FLOAT(I-l)*(MESH(SN, 

C 5)-MESH(FN,5))/FLOAT(STEP)))*'-2+(MESH(N,8)-(MESH(FN,6)+ 

C FLOAT( I - 1 ) * ( ME SH( SN , 6 ) -ME SH( FN , 6 ) ) /FLOAT( STEP ) ) ) **2 ) 

DPS IC=PDVEC ( FN ) +( I - 1 ) *( PDVEC ( SN ) - PDVEC ( FN ) ) / ( STEP ) 

PSIRP = PVEC(FN)+(I-1)*(PVEC(SN)-PVEC(FN))/(STEP) 

DOT = (N0RM11*(MESH(N,7)-(MESH(FN,5)+(I-1)*(MESH(SN,5)- 
C MESH(FN,5))/STEP))+(N0RM12*(MESH(N,8)-(MESH(FN,6)+(I-1)* 

C ( MESH( SN , 6 ) -MESH( FN , 6 ) ) /STEP ) ) ) ) /RMRP 

ENDIF 

CALL HANK CMPLX( RMRP, 0. 0) ,H0RP1,H1RP1) 

IF( ( I . EQ. 1 ) . OR. ( I . EQ. STEP+1 ) ) THEN 

PSI = ( J/4. 0)*(H0RP1*DPSIC - PSIRP*D0K^HlRPl)/2 

ELSE 

PSI = ( J/4. 0)*(H0RP1*DPSIC - PSIRP*D0K^H1RP1) 

ENDIF 

INTEGRAL = INTEGRAL + PSI*DIST/STEP 
CONTINUE 

SUM = SUM + INTEGRAL 
CONTINUE 
T(N,M) = SUM 
CONTINUE 
CONTINUE 

DO 55, I = 1, PERND 

TEMP(I) = (0.0,0. 0) 

DO 50, M = 1, PERND 

TEMP(I) = TEMP(I) + T(I,M)*BCOND(M) 

CONTINUE 

WRITE(2,*) (T(I,MM), MM = 1 , PERND) 

CONTINUE 

DO 60, I = 1, PERND 

WRITE(2,*) I, TEMP(I) 

CONTINUE 

FORMATC IX, 'COLUMN ',13,', ROW ',13,' OUT OF ’,13) 

CLOSE(2) 

RETURN 

END 



SUBROUTINE CNSOLV(T , BCOND , CNVEC , PERND) 

THIS SUBROUTINE CALCULATES THE C VECTOR BY SOLVING: 

-1 

Cn = [I - T] * BOUNDARY CONDITIONS (INCIDENT FIELDS) 



COMPLEX BCOND(100),T(100,100),TEMP(100),CNVEC(100), DETERM 
REAL COND, DMAG 
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INTEGER PERND,I,J,K,L,INORM,NMAX 
C 

INORM = 0 
NMAX =100 
C 

DO 10, 1=1, PERND 

DO 10, J = 1, PERND 
IF(I.EQ. J) THEN 

T(I,J) = 1 - T(I,J) 

ELSE 

T(I,J) = -T(I,J) 

ENDIF 

10 CONTINUE 

WRITE(*,*) 'INVERTING THE [ I - T] MATRIX 
C 

CALL CSMINVC T , NMAX , PERND , DETERM , COND , INORM) 

C 

DMAG = CABS (DETERM) 

WRITE(*,*) ' DETERMINANT = ', DMAG 

WRITE(*,*) ' CONDITION NUMBER = ' , COND 

WRITE(*,*) ' MULTIPLING MATRICES TO FORM THE SCATTERED FIELDS ' 
C 

DO 30, I = 1, PERND 

CNVEC(I) = (0. 0,0. 0) 

DO 30, J = 1, PERND 

CNVEC(I) = CNVEC(I) + T( I , J)*BCOND( J) 

30 CONTINUE 
C 

RETURN 

END 



SUBROUTINE FFLD( CNVEC , PERND , MESH , U , OFFSET , K , NRES , XORG , YORG) 

THIS SUBROUTINE CALCULATES THE FAR FIELDS DUE TO THE OFFSET 
BOUNDARY SCATTERED FIELDS AND THE PERIMETER SCATTERED FIELDS. 
ADDITIONAL GREEN'S FUNCTION INTEGRALS ARE ACCOMPLISHED. 

COMPLEX CNVEC( 100) ,U( 100,100) , J,PSISP( 100) ,TEMP,PSI ,DPSI 
COMPLEX DPSISPC 100) , INTEGRAL, DPS IM( 100, 100) ,PSIM( 100,100) 

REAL MESH( 100, 8), OFFSET, K, DOT, DOTl, PI, ARES, XORG, YORG, DIST, SIGMA 
REAL R,DZ,DL,NORMll,NORM12 
INTEGER PERND , I , M , N , L , NRE S , FN , SN , STEP , 1 1 
C 

WRITEC*,*) ' CALCULATING THE SCATTERED FIELDS' 

J = (0. 0,1.0) 

PI = 4. 0*ATAN(1. 0) 

ARES = 2. 0*PI/FLOAT(NRES) 

STEP = 5 
C 

DO 5, M = 1, PERND 

DO 5, 1=1, PERND 

IF(M. EQ. I) THEN 

PSIM(I,M) = (1.0 + U(I,M))/2.0 
DPSIM(I,M) = (1.0 - U( I, M)) /OFFSET 

ELSE 
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PSIM(I,M) = U(I,M)/2. 0 
DPSIM(I,M) = -U( I, M) /OFFSET 

ENDIF 

5 CONTINUE 
C 

DO 10, I = 1, PERND 

PSISP(I) = (0. 0,0.0) 

DPSISP(I) = (0. 0,0.0) 

DO 10, L = 1, PERND 

PSISP(I) = PSISP(I) + PSIM(I,L)*CNVEC(L) 

DPSISP(I) = DPSISP(I) + DPSIM(I,L)*CNVEC(L) 

10 CONTINUE 

C 

DO 30, 1=0, NRES-1 

WRITE(6,1000) I+l, NRES 
INTEGRAL = (0. 0,0. 0) 

DO 20, N = 1, PERND 
FN = N 

IF(N. EQ. PERND) THEN 
SN = 1 

ELSE 

SN = N + 1 

ENDIF 

R = MESH(SN,5) - MESH(FN,5) 

DZ = MESH(SN,6) - MESH(FN,6) 

DL = SQRT(R**2 + DZ**2) 

NORMll = -DZ/DL 

N0RM12 = R/DL 

DO 15, II = 1, STEP+1 

DIST = SQRT((MESH(FN,5)-MESH(SN,5))**2+(MESH(FN,6)- 
C MESH(SN,6))**2) 

IF(II. EQ. 1) THEN 

PSI = PSISP(FN)+0. 25*(PSISP(SN)-PSISP(FN))/ 

C FLO AT( STEP) 

DPSI = DPSISP(FN)+0. 25*(DPSISP(SN)-DPSISP 
C (FN))/FLOAT(STEP) 

DOT = N0RM11*SIN(I*ARES) + NORM12*COS( I*ARES) 

DOTl = (MESH(FN,5)+0. 25*(MESH(SN,5)-MESH(FN, 
5))/FLOAT(STEP)-XORG)*SIN(I*ARES) + (MESH(FN,6)+ 

0. 25*(MESH(SN,6)-MESH(FN,6))/FLOAT(STEP)- 
YORG)*COS(I*ARES) 

ELSEIFC II. EQ. STEP+1) THEN 

PSI = PSISP(FN)+(FLOAT(STEP)-0.25)*(PSISP(SN)- 
PSISP(FN))/FLOAT(STEP) 

DPSI = DPSISP(FN)+(FLOAT(STEP)-0. 25)*(DPSISP(SN)- 
DPSISPC FN) ) /FLOAT( STEP) 

DOT = N0RM11*SIN(I*ARES) + NORM12*COS(I*ARES) 

DOTl = (MESH(FN,5)+(FLOAT(STEP)-0. 25)*(MESH(SN,5)- 
MESH(FN,5))/FLOAT(STEP)-XORG)*SIN(I*ARES) + (MESH( 
FN,6)+(FLOAT(STEP)-0. 25)*(MESH(SN,6) -MESH(FN,6) )/ 
FLOAT( STEP) -YORG)*COS( I*ARES) 

ELSE 

PSI = PSISP(FN)+FLOAT(II-l)*(PSISP(SN)-PSISP(FN))/ 
C FLOAT(STEP) 

DPSI = DPSISP(FN)+FL0AT(II-1)*(DPSISP(SN)-DPSISP 
C (FN))/FLOAT(STEP) 
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DOT = N0RM11*SIN(I*ARES) + N0RM12*C0S( I*ARES) 

DOTl = (MESH(FN,5)+FL0AT(II-1)*(MESH(SN,5)-MESH(FN, 
C 5))/FL0AT(STEP)-X0RG)*SIN(I*ARES) + (MESH(FN,6)+ 

C FL0AT( I I - 1 )*( MESH( SN , 6 ) -MESH( FN , 6 ) ) /FL0AT( STEP ) - 

C Y0RG)*C0S(I*ARES) 

ENDIF 

IF((II.EQ. 1). OR. (II. EQ. STEP+1)) THEN 

TEMP=( J*DPSI+DOT*PSI)*(EXP( J*D0T1) )*DIST/(2. 0* 
C STEP) 

ELSE 

TEMP=( J*DPS I+DO'P'^PS I )*( EXP( J*D0T1 ) )*DI ST/ 

C (STEP) 

ENDIF 

INTEGRAL = INTEGRAL + TEMP 
15 CONTINUE 

20 CONTINUE 

SIGMA = ( (CABS( INTEGRAL) )**2. 0)/(4. 0*K) 

WRITE(50,1010) I+l, (I*ARES*180. 0/PI), SIGMA 
30 CONTINUE 
C 

WRITE(6,*) ' ’ 

WRITE(6,*) ' «< FIELD PATTERN STORED IN FFPAT.DAT »> ' 

WRITE(6,*) ' ’ 

C 

1000 F0RMAT( IX, 'INTEGRAL ',13,', OUT OF ’,13,' COMPLETED’) 

1010 F0RMAT(1X,I3,2X,F6. 2,2X,E14. 8) 

C 

RETURN 

END 

C 

C 

SUBROUTINE HAN1(Z,H0,H1) 



Computing Hankel Functions for n=0,l with 
Complex Argument, Z. Direct Power Series Method for 
CABS(Z) . LE. 5 and Hankel's Asymptotic Formula for 
CABS(Z) .GT. 5. Written 11/6/87 by M. A. Morgan 

INTEGER M,M2 

REAL C(34) ,DM,F(34),G0,P(34),Pi,P2 

COMPLEX Z,Z2,Z3,Z4,J0, J1,Y0,Y1,AM,CL,P0,P1,Q0,Q1 

COMPLEX E0,E1,X0,X1,H0,H1, j 

PI=3. 1415927 

P2=2. 0/PI 

j=(0. ,1. ) 

IF(CABS(Z). LE. 5. 0) THEN 

Direct Power Series Method 

G0= 1. 78072 
Z2=0. 5*Z 
CL=CL0G(G0*Z2) 



Computing F(m) = m ! and P(m) = 1 + 1/2 + 1/3 + ....+ 1/m 
F(l)=l. 0 



130 



P(l) = l. 0 
DO 11 M=2,34 

F(M)=M*F(M-1) 

P(M)=P(M-1)+1. 0/M 
11 CONTINUE 

C 

C Computing Power Series Coefficients 
C 

DM=-1. 0 
DO 22 M=l,34 

C(M)=DM/(F(M)*F(M)) 

DM=-DM 

22 CONTINUE 

C 

C Computing JO and J1 
C 

J0=(1. ,0. ) 

J1=(0. ,0. ) 

M=0 

33 M=M+1 

M2=2*M 

AM=C(M)*(Z2**M2) 

J0=J0+AM 

J1=J1-M*AM 

IF((CABS(AM).GT. 1.0E-10).AND. (M. LT. 34)) GO TO 33 
J1=J1/Z2 
C 

C Computing YO and Y1 
C 

M=0 

Y0=CL*J0 

Yl=Z2*CL*Jl-0. 5*J0 
44 M=M+1 

M2=2*M 

AM=C ( M ) *P ( M ) * ( Z 2**M 2 ) 

Y0=Y0-AM 

Y1=Y1+M*AM 

IF((CABS(AM). GT. 1. OE-10). AND. (M. LT. 34)) GO TO 44 

Y0=P2*Y0 

Y1=P2*Y1/Z2 

H0=J0-j*Y0 

Hl=Jl-j*Yl 

RETURN 

ELSE 

C 

C Hankel' Asymptotic Formula (Abram. & Stegun p. 364 
C 

Z2=Z*Z 

Z3=Z*Z2 

Z4=Z*Z3 

P0=1. 0-. 0703125/Z2+. 1121521/Z4 
Q0=-. 125/Z+. 0732422/Z3 
Pl=l. 0+. 1171875/Z2-. 1441956/Z4 
Ql=. 375/Z-. 10253906/Z3 
X0=(Z-. 25*PI) 

X1=(Z-. 75*PI) 
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E0=CEXP(-j*X0) 
E1=CEXP( -j*xi) 
AM=CSQRT(P2/Z) 
H0=AM*(P0-j*Q0)*E0 
Hl=AM*(Pl-j*Ql)*El 

END IF 
C 

RETURN 

END 



SUBROUTINE CSMINV( A, NDIM,N, DETERM, COND, INORM) 

INORM - FLAG TO NORMALIZE COLUMNS AND ROWS OF MATRIX A 



MATRIX 


NORMALIZATION BY M. A. 


MORGAN 




APRIL : 


24,1978 






A 


- MATRIX TO INVERT 




- INPUT/ OUTPUT 


NDIM 


- 




- INPUT 


N 


- 




- INPUT 


DETERM 


- DETERMINATE OF A 




- OUTPUT 


COND 


- CONDITION NUMBER OF 


A 


- OUTPUT 


INORM 


- INTEGER NORMALIZATION FLAG 


- INPUT 



COMPLEX A( 100, 100) ,PIVOT(100) ,AMAX,T, SWAP, DETERM, U 
INTEGER I,J,K,L,IPIVOT(100) ,INDEX(100,2) ,IROW,ICOLUM,Ll,JROW 
INTEGER JCOLUM,N, INORM 

REAL TEMP,ALPHA( 100) ,COL( 100) ,ROW( 100) ,AJK, SUMAXA, SUMROW, SUMAXI 
IF(NDIM. GT. 100) THEN 

WRITE(*,*) ' ERROR IN INVERTION CALL... DIMENSION > 100 ' 
STOP 

END IF 

IF(N. GT. NDIM) THEN 

WRITE(*,*) ' ERROR IN INVERTION CALL... N > MAX DIM. ' 

STOP 

ENDIF 

IF( INORM. NE. 1) GO TO 7 
DO' 3 K = 1, N 

COL(K) =0.0 
DO 1 J = 1,N 

AJK = CABS(A(J,K)) 

IF(AJK. GT. COL(K)) COL(K) = AJK 
CONTINUE 
DO 2 J = 1, N 
A(J,K) = A(J,K)/COL(K) 

CONTINUE 
ROW NORMALIZING 
DO 6 J = 1, N 

ROW(J) = 0.0 
DO 4 K = 1, N 

AJK = CABS(A(J,K)) 

IF(AJK. GT. ROW(J)) ROW(J) = AJK 
4 CONTINUE 
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DO 5 K = 1, N 

5 A(J,K) = A(J,K)/ROW(J) 

6 CONTINUE 

7 CONTINUE 

DETERM = CMPLX(1. 0,0. 0) 

SUMAXA = 0. 0 
DO 20 J = 1, N 

ALPHA(J) = 0.0 
SUMROW = 0. 0 
DO 10 I = 1, N 

ALPHA(J) = ALPHA(J) + A(J,I)* CONJG(A( J, I) ) 

10 SUMROW = SUMROW + CABS(A(J,I)) 

ALPHA(J) = SQRT(ALPHA(J)) 

IF( SUMROW. GT. SUMAXA) SUMAXA = SUMROW 
20 IPIVOT(J) = 0 

DO 600 I = 1, N 
C 
C 

AMAX = CMPLXCO. 0,0. 0) 

DO 105 J = 1, N 

IF (IPIVOTC J)-l) 60, 105, 60 
60 DO 100 K = 1, N 

IF (IPIVOT(K)-l) 80, 100, 740 

80 TEMP = AMAX*CONJG(AMAX) - A( J,K)*CONJG(A( J,K) ) 

IF(TEMP)85,85,100 

85 IROW = J 

ICOLUM = K 
AMAX = A(J,K) 

100 CONTINUE 

105 CONTINUE 

IPIVOT( ICOLUM) = IPIVOTC ICOLUM) + 1 
C 
C 

IF (IROW- ICOLUM) 140, 260, 140 
140 DETERM = -DETERM 

DO 200 L = 1, N 

SWAP = A(IROW,L) 

A(IROW,L) = A( ICOLUM, L) 

200 A( ICOLUM, L) = SWAP 

SWAP = ALPHA(IROW) 

ALPHA(IROW) = ALPHAC ICOLUM) 

ALPHAC ICOLUM) = SWAP 
260 INDEX(I,1) = IROW 

INDEX(I,2) = ICOLUM 
PIVOT(I) = A( ICOLUM, ICOLUM) 

U = PIVOT(I) 

DETERM = DETERM*U 

DETERM = DETERM/ ALPHAC ICOLUM) 

TEMP = PIVOT(I)*CONJG(PIVOT(I)) 

IF(TEMP)330,720,330 

C 

C 

330 A( ICOLUM, ICOLUM) = CMPLX( 1. 0 , 0. 0) 

DO 350 L = 1, N 
U = PIVOT(I) 

350 A( ICOLUM, L) = A( ICOLUM, L)/U 
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c 

c 

380 

400 



450 

550 

600 

C 

C 

620 



630 



705 

710 



900 

910 



950 

955 

720 

730 

740 

C 

C 



DO 550 LI = 1, N 

IF(Ll-ICOLUM) 400, 550, 400 
T = A(L1,IC0LUM) 

A(L1,IC0LUM) = CMPLX(0. 0,0. 0) 
DO 450 L = 1, N 

U = A(ICOLUM,L) 

A(L1,L) = A(L1,L) - U*T 

CONTINUE 

CONTINUE 



DO 710 I = 1, N 

L = N + 1 - I 

IF (INDEX(L,1) - INDEX(L,2)) 630, 710, 630 
JROW = INDEX(L,1) 

JCOLUM = INDEX(L,2) 

DO 705 K = 1, N 

SWAP = A(K,JROW) 

A(K,JROW) = A(K, JCOLUM) 

A(K, JCOLUM) = SWAP 
CONTINUE 

CONTINUE 
SUMAXI = 0.0 
DO 910 I = 1, N 
SUMROW = 0. 0 
DO 900 J = 1, N 
SUMROW = SUMROW + CABS(A(I,J)) 
IF(SUMROW.GT. SUMAXI) SUMAXI = SUMROW 
CONTINUE 

COND = SUMAXA*SUMAXI 
IFCINORM. NE. 1) GO TO 955 
DO 950 K = 1, N 

DO 950 J = 1, N 

A(J,K) = A(J,K)/(ROW(K)*COL(J)) 

CONTINUE 

RETURN 

WRITE(*,730) 

FORMATC ' MATRIX IS SINGULAR’) 

RETURN 

END 
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APPENDIX F. DIELECTRIC CYLINDER SCATTERING PROGRAM 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 



PLANEWAVE SCATTERING BY A DIELECTRIC CYLINDER 
E - WAVE (TM CASE) 

H - WAVE (TE CASE) 



CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 



c 



c 



c 



c 



c 



c 



c 



c 



COMPLEX GAMMA(0: 200,2) ,SIGMA(2) ,ER,MU,KR,YR, ZR 
COMPLEX*16 JA(0; 200) ,DJA(0: 200) ,KRRA 
COMPLEX*16 J(0: 200) ,DJ(0: 200) 

REAL*8 Y(0: 200) ,DY( 0: 200) ,YA(0: 200) ,DYA( 0: 200) , JB( 0: 200) 
REAL*8 DJB(0: 200) ,RA, KORA 
REAL K0,PI,SIGMAN(2),A,B 
INTEGER I, I I, MODE, ARES, N 



OPEN(3,FILE='C: MSFORT DECTEM.DAT’) 
PI = 3. 1415927 



WRITE(*,*) 'PERMITTIVITY FORMAT IS " a + jb, " ' 

WRITE(*,'') 'Enter Dielectric Constant (REAL PART, a) ' 
READ(*,*) A 

WRITE(*,*) 'Enter Dielectric Constant (IMAGINARY PART, b) ' 
READ(*,*) B 
ER = CMPLX(A,B) 



WRITE(*,*) 'PERMEABILITY FORMAT IS " a + jb, " ' 

WRITE(*,") 'Enter Permeability Constant (REAL PART, a) ' 
READ(*,*) A 

WRITE(*,*) 'Enter Permeability Constant (IMAGINARY PART, b) ' 
READ(*,*) B 
MU = CMPLX(A,B) 



WRITE(*,*) 'Enter the wave number (ko) ' 
READ(*,*) KO 



WRITE(*,*) 'Enter Cylinder Radius (IN WAVENUMBER UNITS) 
WRITE(*,*) 'WARNING: Do Not Enter Zero ! ' 

READ(*,*) KORA 

KR = CSQRT(MU*ER) 

YR = CSQRT(ER/MU) 

ZR = CSQRT(MU/ER) 

RA = KORA/KO 
KRRA = KR*RA 

WRITE(*,*) 'Enter No. of Modes: ' 

READ(*,*) MODE 

WRITE(*,*) 'Enter the angular resolution ' 

READ(*,*) ARES 
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o o o n n 



c 



c 

c 

c 



10 

c 

c 

c 



20 



40 

30 

C 

C 

110 



120 

1000 

C 

C 



WRITE(3,*) 'Cylinder Scattering vs. angle' 
WRITE( 3 , 1 1 0 ) ER , MU , RA , MODE , KO , KRRA 
WRITE (3,*) ' ' 



CALL BES( MODE, KORA, JB,Y,DJB,DY) 

CALL DCBJNS ( DCMPLX( KORA, 0. OD+00) , MODE, J,DJ) 

CALL DCBJNS ( (KRRA*K0) ,MODE, JA,DJA) 

WRITE(*,*) ' RETURNED FROM FINAL BESSEL CALL' 

CALCULATING GAMMAs 

DO 10, N = 0, MODE 

GAMMA(N, 1) = ( JA(N)*DJ(N)-YR*DJA(N)*J(N))/(JA(N)*CMPLX(DJ(N) 
C ,-DY(N)) - YR*DJA(N)*CMPLX(J(N) ,-Y(N))) 

GAMMA(N,2) = ( JA(N)*DJ(N) -ZR*DJA(N)*J(N) )/( JA(N)*CMPLX(DJ(N) 
C ,-DY(N)) - ZR*DJA(N)*CMPLX(J(N),-Y(N))) 

WRITE( * , 1000 ) N , GAMMA( N , 1 ) , GAMMA( N , 2) 

CONTINUE 

CALCULATING SIGMAs 

DO 30, I = 180, -180, -ARES 
DO 40, II = 1, 2 

SIGMA(II) = (0. 0,0. 0) 

DO 20, N = 1, MODE 

SIGMA(II) = SIGMA(II)+2. 0*GAMMA(N,II)*COS(N*PI*I/180. 0) 
CONTINUE 

SIGMA(II) = SIGMA(II) + GAMMA(0,II) 

SIGMAN(II) = ((4. 0/K0)*(CABS(SIGMA(II)))**2) 

CONTINUE 

WRITE(3,120) I, SIGMAN(l), SIGMAN(2) 

CONTINUE 



FORMATdX, 'Er = ' ,F6. 4, 1X,F6. 4,/, 



c' 


Mu = 


' ,F6.4,1X,F6.4,/, 


c' 


RADIUS ( meters) = 


' ,F8.5,/, 


c' 


MAX MODE 


’ ,13,/, 


c' 


KO 


' ,F8.5,/, 


c' 


KrRa 


' ,F8. 5,1X,F8. 5) 



F0RMAT(1X,I4,2(3X,E14.4)) 
F0RMAT(1X,I3,1X,2(E12.4,1X,E12. 4,4X)) 



STOP 

END 
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APPENDIX G. SOFTWARE SOURCES 



1. DISSPLA 

Integrated Software Systems Corporation 
10505 Sorrento Valley Road 
San Diego, CA 92121 

2. CURVE-DIGITIZER 
West Coast Consultants 

4202 Genesee Avenue, Suite 309 
San Diego, CA 92117 

3. Microsoft FORTRAN 
16011 NE 36th Wav 
BOX 97017 
Redmond, WA 98073 

4. Microwav NDP FORTRAN 
POB 79 

Kingston. MA 02364 

5. Prof MA. Morgan. Code 62Mw 
Naval Postgraduate School 
Monterey, CA 93943 

6. LT T.B. Welch 111 
c 0 T.B. Welch Jr. 

1318 Walthour Road 
Savannah, GA 31410 
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